Self-evolving atomistic kinetic Monte Carlo: fundamentals and applications

The evolution of defects in materials is a major topic in physics, chemistry and materials science [1]. It is intrinsically a challenging multiscale problem in time and length. As a valuable tool, computer simulation has been widely employed to study the equilibrium and transport properties of defects [2]. For instance, molecular dynamics (MD) is extensively used in investigating atomistic processes and thermodynamic properties of materials [3]. However, due to the required time step, the accessible time scale in MD is generally limited to nanoseconds or less. Consequently, only short-term defect evolution can be directly simulated using MD. Long-term defect evolution, which is vital in determining many materials properties, such as irradiation induced segregation or void swelling [4], is beyond the scope of MD and requires other techniques. Several approaches have been suggested to extend the time scale and bridge the gap between atomistic simulations and experiments. Accelerated dynamics [5] and kinetic Monte Carlo (KMC) [6] are two such methods.

Accelerated dynamics approaches employ a variety of techniques to speed up the atomic transition rate while maintaining the framework of MD. For instance, parallel replica dynamics [7] explores the phase space using multiple trajectories at the same time. Temperature accelerated dynamics (TAD) [8] increases the simulation temperature to boost the possibility of finding escape pathways and advance time at a rate appropriate to overcoming the lowest barrier at the low target temperature. Hyperdynamics [9] uses a bias potential to raise the potential minima while keeping the saddle points (SPs) unchanged. Similar to hyperdynamics, metadynamics [10] introduces Gaussian functions to the potential energy surface to force the atoms to explore other regions of phase space. An advantage of these methods (except for metadynamics) is that the true dynamics is theoretically followed. However, many difficulties remain for practical applications. For instance, a general method of specifying the bias potential for hyperdynamics is currently unavailable [5]. Moreover, the scaling of these algorithms generally limits system size, which prevents certain processes and phenomena being investigated.

As an alternative to accelerated dynamics, KMC is a stochastic method and has been used to study defect evolution for decades [6]. Conventional KMC includes object KMC (OKMC) [11] and atomistic KMC (AKMC) [12]. In OKMC, defects are considered as abstract objects and no atomistic details are included. This leads to a significant improvement in computational speed. However, it also causes two major limitations of OKMC: (i) the saddle points cannot depend on the defect configuration, and (ii) the defect–defect interaction cannot be properly described. Most current AKMC models employ an on-lattice approximation, which limits their ability to describe a system in which defects are not exactly on-lattice (e.g. interstitials) or systems which undergo significant deformation. Therefore, AKMC is generally limited to studies of vacancies [12]. For both OKMC and AKMC, defining a list of possible mechanisms a priori is required. This is a severe limitation of the conventional KMC methods, since atomic motion may not be intuitive and can be extremely difficult to guess or determine in advance [5].

To overcome these limitations, 'on-the-fly' approaches based on an off-lattice system have been developed, such as adaptive KMC [13, 14], the kinetic activation–relaxation technique (KART) [15], autonomous basin climbing (ABC) [16] and self-evolving atomistic KMC (SEAKMC) [17]. The common feature of these methods is that saddle points are found 'on the fly' during the simulation in an off-lattice system. The difference resides in the ways of searching for saddle points and other details. Adaptive KMC uses the dimer method [18]; kinetic ART employs the Lanczos method [19], ABC uses autonomous basin climbing [20, 21] with a nudged elastic band (NEB) [22, 23] and SEAKMC employs dimer search in combination with selective active volumes (AVs) [17]. These approaches provide significant improvement over the conventional KMC. However, each method has its own limitations. For instance, ABC only uses the minimum energy saddle point in KMC for each defect, which limits the system evolution with a bias towards a ground state and prohibits the high energy fluctuations in states. In addition, several common issues remain. With the exception of SEAKMC, which has been applied to cascade annealing simulations with ∼128 000 atoms, the application of all these methods has been limited to small system sizes and very few barriers, although the capability to study long-term defect evolution with atomistic fidelity in a system with multimillion atoms is desirable in many areas. In addition, the error estimations in both saddle points and simulation time have not been properly carried out.

In this paper, we present the further development of SEAKMC as well as its applications. In particular, the method of determining the AVs and the effect of AV size on saddle point energies are elucidated. A multiple-step procedure involving the use of a reduced AV size initially, followed by the desired AV, is proposed to increase the efficiency of the method. For saddle point searches, the observation that the dimer method has a higher probability of finding the saddle point with lowest energies is investigated and the saddle point distributions for various defects are established. In order to estimate the fidelity of dynamics and accuracy of the simulation time, a general criterion is proposed and applied to two representative scenarios. Applications of SEAKMC for investigating the diffusion of self-interstitial atoms (SIA) and vacancies in bcc iron are presented and compared with MD simulations, which demonstrates that the kinetic Monte Carlo method can be used to investigate dynamic processes and generate results that formerly could be obtained only through MD. Furthermore, the correlation factors for vacancy and SIA diffusion have been calculated using SEAKMC, which proves that SEAKMC possesses some unique predictive capabilities. The strengths and limitations of SEAKMC are discussed to present a comprehensive picture of the SEAKMC method in both technically important details and its potential broad applications.

The framework of SEAKMC consists of several components: selection of AVs, saddle point searching, kinetic Monte Carlo and static relaxations [17]. A flowchart for the SEAKMC method is given in figure 1. The SEAKMC first identifies an AV for each defect. Then saddle points are searched for within each AV, followed by a KMC step, which randomly chooses a particular saddle point based on probabilities and advances the simulation time. Depending on the defect nature and arrangements of AVs, either a local (within the AV) or a global (the whole system) static relaxation is carried out to move the system from the saddle point to another local minimum. The procedure repeats itself until the desired statistical results can be extracted or the desired simulation time has been reached. The details of each component mentioned above are given in this section.

Zoom In Zoom Out Reset image size

2.1. Active volume

The concept of AV was introduced [17] based on the fact that most defects and their saddle points are localized. The atoms far away from the defects or dynamic processes are essentially irrelevant and have little effect on these defects or processes. Correspondingly, the perturbation caused by the existence of defects on the atoms far away is also negligible and the properties of these atoms are extremely close to the case in which no defect is present. Therefore, an AV (localized region) containing one or more defects can be defined according to several criteria. In contrast to a previous method [14], the method focuses on this part of the system not only for the initial displacement but also the following dimer searches. All the atoms outside the AVs are frozen during the dimer searches. The effect of internal defects should be negligible at the AV boundary, and atoms outside the AVs should have no significant effects on the internal defects or processes. An alternative way to express this is that the AVs should be sufficiently large that the properties of internal defects can be reasonably described by considering only the atoms inside AVs. To practically establish the AV boundaries, several possible criteria could be used, such as the maximum allowed deviation in energy, strain or stress relative to the perfect bulk value. When the AV considered is under the influence of a long-range strain field, such as from a dislocation, the criteria adopted need to account for the atomic strain.

Figure 2 shows the energy difference from the bulk value for an atom as a function of the distance between the atom and a defect. The two types of defect considered here are a self-interstitial atom (SIA) in the dumbbell configuration and a vacancy. It is clear that the energy difference gets smaller as distance increases and the vacancy is more localized than the dumbbell since the energy difference approaches zero at a shorter distance for the vacancy (4 versus 6 lattice constants). This means that the AV size could be different for different defects even when the same criterion is used. For example, when an energy deviation of 0.001 eV is chosen as the criterion, the AV configurations for dumbbell and vacancy are shown in figure 3. Combining figures 2 and 3, it can be seen that the AV can be potentially defined as a sphere with a radius of a few lattice constants. To test this, the dependence of saddle point energy on AV size was investigated and the results are shown in figure 4. The saddle point energies could significantly deviate from the correct value when the AV size is small. For instance, the migration energy is 0.377 eV when the radius of the AV is 2.63 lattice constants while the bulk value is 0.304 eV using the Ackland-04 interatomic potential [24]. Nevertheless, the saddle point energies approach the correct value when the AV size is larger. Moreover, consistent with the above analysis, the value for a vacancy approaches the bulk value at a shorter distance than it does for the dumbbell. To briefly summarize, analyzing the energy difference from the bulk value of the atoms, a certain criterion can be chosen to define the AV size. By establishing the dependence of transition energies on AV size, a proper AV size can be determined, which provides information back to the criterion selection process. Therefore, these two steps are complementary and are both essential in defining AVs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The efficiency of SEAKMC is mainly a result of the calculation focusing on AVs. Therefore, the size of AVs has a significant effect on the efficiency of the method. For instance, choosing a smaller deviation of energy from the bulk value will lead to a larger AV size. This results in a higher accuracy but less acceleration. In practice, the optimum criterion is usually chosen as a compromise between accuracy and computational expense. The shape of AVs discussed so far has been spherical due to symmetry of the internal defects. However, in general, AVs could have any shape, depending on the defect properties of the crystal lattice in the vicinity of the AVs. For instance, the AV could be cylindrical for a dislocation line and ellipsoidal for small glissile interstitial clusters. The optimization of the AV shape and size provide additional enhancement of the SEAKMC. The efficiency of SEAKMC is estimated to increase as system size gets larger, in particular for the case where defects are well separated. This is because the vast majority of the calculations involve only a small part of the system volume, i.e. a few AVs. Therefore, as few as a few hundred or thousand atoms may be computationally active in a multimillion atom system. The efficiency will decrease if defects are close to each other or if long-range defect interactions need to be considered, since larger AVs are required. However, the computational cost is still considerably lower than using the whole system.

In addition, the criterion for selecting AVs or the AV size can evolve during the simulations to provide further acceleration. The most time-consuming part of SEAKMC is the saddle point searching. Each search generally involves hundreds to thousands of force calculations. We suggest a multi-step procedure which begins with a relatively small AV that is used to obtain an initial prediction of the transition states. Then, the AV size is gradually increased to obtain more accurate values. The efficiency of SEAKMC in this scheme is enhanced significantly, not only because the forces need to be calculated for fewer atoms, but also because the number of interactions required to obtain the saddle points is reduced since fewer degrees of freedom are involved. The saddle point configurations obtained in the reduced AV are usually close to the actual ones and the further corrections to the saddle point energy and configuration of each saddle point are made then in the larger AV with minimum computational cost. In practice, the AV size is increased until the saddle point energies converge to the desired accuracy. The efficiency of this multi-step procedure is found to be higher than using a large AV directly while maintaining the same accuracy. For example, using a 10 × 10 × 10 system with 2000 atoms for bcc iron, the speedup using this multi-step procedure for a single vacancy and SIA is found to be 115 and 37 times faster, respectively, than using the whole system for dimer search. Furthermore, the boost is expected to be more significant as the system size increases because the AV is a small fraction of the volume. It is noted that starting from a too small AV and with fewer degrees of freedom may prevent some saddle points being found, although there is a possibility that such saddle points can be recovered during the search in the larger AVs. This procedure is absolutely safe when the defect or process is localized, such as diffusion of point defects, but may cause some additional inaccuracy if the saddle points involve coordinated movements of a large number of atoms. Therefore, depending on the nature of the problem, caution is needed to determine whether this further acceleration procedure should be employed.

For the examples presented here, i.e. bulk self-interstitial atom diffusion, we have used a simplified and faster procedure where the radius of a spherical AV was estimated from the data presented in figures 2 and 3. In this particular case the local environment around the vacancy is always the same and a spherical AV of the same radius can be applied after each vacancy jump. However, if several defects are included in the simulated system, their AVs will evolve according to the interactions between them. For instance, if defects are close to each other and could interact, they will be included in a single AV so that the collaborative movement of the defects can be included during the saddle point searches, which may lead to irregular shapes of AVs if energy difference is chosen as the criterion. Once the system moves to the corresponding saddle point, the relaxation will be carried out to include the possible defect interactions.

2.2. Saddle point searching

The saddle point searches can be carried out using several techniques, such as dimer, Lanczos and ABC with NEB [25]. Dimer and Lanczos search for saddle points directly, without knowledge of the final state a priori. ABC finds the final states then uses NEB to determine the saddle points. It should be noted that the dimer, Lanczos and NEB methods are based on the harmonic approximation of transition state theory (HTST) [18, 19]. Therefore, all current 'on-the-fly' KMCs are similarly dependent on HTST. Since ABC is not suitable for finding multiple saddle points simultaneously, especially for saddle points with the same energy barrier, it is not fully consistent with KMC since KMC requires a full catalog of all the possible saddle points in order to accurately calculate time. Neglecting certain saddle points introduces error in the KMC time. In this study, the dimer method is chosen to illustrate the concept. Comparisons with the Lanczos method are made wherever possible.

We have implemented two different versions of the dimer method, the improved dimer [18, 26] and superlinearly converging dimer [27]. The details of these two approaches can be found in previous publications [18, 26, 27]. We briefly mention the difference here. The improved dimer method uses the conjugate gradient (CG) method for dimer rotation and translation while the superlinear dimer method employs the limited memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) technique. For applications in this study, no qualitative difference was observed when using these two methods. For instance, if the same system setup is used, both methods predict the same energy barriers. Furthermore, no significant improvement of convergence was found using the superlinearly converging dimer method, in contrast to a previous study [27]. The improved dimer and superlinear dimer methods require similar numbers of iterations for the system to reach the saddle points when studying interstitials and vacancies in Fe. Since each saddle point search is independent of the others, the SEAKMC code can be parallelized using the hybrid MPI/OpenMP framework [28].

The saddle point distribution calculated for the SIA and vacancy using two different interatomic potentials with the dimer method is given in figure 5. The SIA results include both the [110] dumbbell and [111] crowdion configurations. The saddle point with the lowest energy have the highest possibility of being found, which is similar to results reported in previous studies [13]. To determine whether this is a general feature of the method, the saddle point distribution for small SIA clusters (di-interstitial (I2), tri-interstitial (I3), tetra-interstitial (I4) and penta-interstitial (I5)) were also calculated as shown in figure 6. The di-interstitial exhibits this behavior while the larger interstitial clusters (I3–I5) do not. Further examination suggests that such behavior is associated with the symmetry of the defect. For instance, the SIA, I2 and vacancy all have relative high symmetry. The saddle point distribution of interstitial clusters obtained in this study is consistent with results calculated using Lanczos method [29]. The distribution function for interstitial clusters is complicated and can be very different for different initial configurations. An SIA cluster of a given size can exist in many different configurations and the results in figure 6 were obtained using two different initial configurations for each cluster size. Because of the number of configurations, the number of distinct saddle points strongly increases with size. This will create a challenge when trying to calculate the simulation time; this issue is discussed further in section 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given how AVs are employed in SEAKMC, recycling of saddle points is relatively straightforward. If the boundary of an AV overlaps another AV that is chosen by KMC to undergo a particular event, the saddle points of this AV should be researched. For the rest of the AVs, saddle points can be directly recycled, since no configuration change occurs in these AVs. Using this scheme, only one or a few AVs are active at any one time in most cases for a multiple-AV system, which increases the efficiency of the inner loop of the SEAKMC method which is indicated in figure 1. Compared with the recycling method used in a previous study, this removes the possibility of using an unphysical initial configuration. The previous study also involved a fixed distance threshold Δr (Δr = 0.2a in [13]) to determine the saddle point configuration combination for recycling. This may cause an artificial strain in the system and lead to a bias in the defect evolution. In contrast, the recycling scheme in SEAKMC is always based on the current state of a relaxed system.

2.3. Kinetic Monte Carlo and static relaxation

Limitations and potential errors of SEAKMC are discussed in this section. In particular, we focus on the fidelity of the dynamics and the accuracy of the simulated physical time. In general, these two issues depend on how complete is the list of distinctive saddle points obtained for each defect configuration. For practical applications, the probability of finding the lowest saddle point and the saddle point distribution are unknown. The configuration of the defect can be extremely complex so that the number of possible saddle points is large, such as for the interstitial clusters mentioned above (see figure 6). In this case, it is extremely difficult to simulate diffusion processes accurately and to determine the simulation time. The fundamental difficulty comes from the limitation of our knowledge about the full list of saddle points. Sampling all the possible saddle points for these complicated defect configurations would require excessive computing; e.g., note the increasing number of saddle points for I4 and I5 in figure 7. The number of unique saddle points continues to increase without saturation for up to 5000 dimer searches, which indicates that it is not practical to require sampling all saddle points for these clusters. This problem is common for all other 'on-the-fly' KMC approaches such as adaptive KMC or kinetic ART. A reasonable and self-consistent scheme for determining when the dimer search is adequately saturated in SEAKMC was formulated based on the incremental change in frequency (Δf) of the events over a given interval (Δn) of dimer searches. This is expressed in the following ratio:

$$\begin{eqnarray} &&\delta ={f}_{\mathrm{n}}^{\ast }/{f}_{\mathrm{ total}} \end{eqnarray} \tag{ 1 }$$

where δ is the ratio of the change in frequency obtained in a specified number of dimer searches, ${f}_{n}^{\ast }=\Delta f/\Delta n$, to the total frequency, f_total, obtained from all the dimer searches. For example, δ can be based on frequency sampled between 4000 and 4100 dimer searches divided by the total frequency. It is noted that this stopping criterion is based on the sampled frequencies, not the number of saddle points. The frequencies are calculated from the migration energy barriers and prefactors, which weights the saddle points differently. The results of applying this formula for various defects are shown in figure 8 for the increment Δn = 5. It is seen that the change in frequency becomes reasonably small for all the cases investigated. However, this scheme cannot completely resolve the fundamental limitation mentioned previously since the error in simulation time cannot be estimated without knowing all the important saddle points. In these situations, we focus on the fidelity of the dynamics. For interstitial clusters using this scheme, we found the dynamics is very close to that in MD. The detailed analysis will be presented elsewhere. The saddle point distribution for interstitial clusters is shown in figure 6, in which A and B stand for two different initial cluster configurations. The previous observation that the lowest saddle point has the highest probability of being found is not true in general. Furthermore, more than 100 different configurations were investigated to confirm this. However, it is seen that many low energy states exist and the energy difference between states is relatively small. Therefore, even if there is a possibility that the saddle point with lowest energy may not be found, the dynamics is observed to be reasonably accurate since a saddle point with very low energy will be chosen in the KMC step with a higher probability. As a result, the dynamics from statistical sampling is not very far from the true trajectory and the probability to realize high energy states still exists. In addition, since the saddle point energies are very small for glissile interstitial clusters, MD is more suitable for simulating their diffusion. In contrast, if all the saddle point energies are high and the temperature is low enough, the SEAKMC will be far more efficient. Thus, the choice of the method depends on the nature of the problem.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If the saddle point distribution is relatively simple and the probability of finding the saddle point with the lowest saddle point energy is higher than finding the others, both the dynamics and simulation time can be correctly described with a relatively low number of saddle point searches. This is usually the case for defects with high symmetry, such as a single vacancy or SIA. Let us consider the accuracy in the simulated defect dynamics and the correlation factor of defect jumps as a measure of the overall fidelity. Figure 9 shows the average cosine value of the angle between two consecutive jumps, from which the defect correlation factor can be extracted using the following equation:

$$\begin{eqnarray} &&f=[1+\langle \cos \theta \rangle ]/[1-\langle \cos \theta \rangle ] \end{eqnarray} \tag{ 2 }$$

where f is the correlation factor and θ is the angle between consecutive jumps. 〈cosθ〉 is the mean cosine averaged over a large number of jumps. We have simulated two different cases of vacancy diffusion in bcc Fe. In the one we have permitted the dimer search to identify all eight distinctive saddle points and their barriers and then applied the usual KMC to choose the particular jump direction. In the other we have halted the dimer search after it identified only one of the saddle points and advanced the process accordingly. The correlation factor in both cases was found to be the same. Figure 9 shows that after ∼2500 jumps 〈cosθ〉 saturates at a value close to zero, which corresponds to a correlation factor of unity, an indicator of random-walk diffusion. A and B indicate two independent simulations of diffusion processes. This confirms that, in the case when there is only one type of saddle point, the first nearest neighbor jump, use of only one random saddle point reproduces the random-walk feature of the diffusion and is equivalent to a full KMC treatment of all the equivalent saddle points. This is because the dimer method with a random search is observed to find equivalent saddle points with equal probabilities. The difference will be in the physical time simulated, and in the case of a single saddle point simulation the time must be corrected for the neglected probability of seven other jump directions. Furthermore, it is found that the same correlation factor is obtained whether one randomly chooses a total of eight saddle points, which may include multiple selections of any of the eight distinct possibilities, or if the eight actual saddle points are selected. Therefore, if the probability of finding the lowest saddle point is known, the dynamics and simulation time can be accurately determined by setting the product of this probability and number of saddle point searches to be the same as the number of required saddle points. For example, the probability of finding the saddle point corresponding to dumbbell diffusion in bcc iron is about 65% (figure 5). Therefore, by setting the number of saddle point searches to be 12, there is a high probability that eight such saddle points will be found, which is the number of first neighbor jumps. This will lead to relatively correct dynamics and simulation time because of the particular distribution of saddle points and their energies. The simulation efficiency is significantly improved by randomly choosing saddle points, compared with guaranteeing that all the required distinct configurations are found. In contrast, the number of unique saddle points found by the dimer method for small interstitial clusters as a function of the number of searches is shown in figure 7. Note that the number of unique saddle points is significantly less than the number of dimer searches. This is because the same saddle point has been found many times during the simulation due to the stochastic nature of the dimer search, and a significant portion of searches lead to the same saddle points. The fact that the probability of finding the same saddle point is not trivial is also a limitation of the approach proposed previously [13]. As described above, it is possible to obtain a sufficiently accurate sampling of the saddle points for single vacancies and interstitials and, since many diffusion studies involve only point defects in pure metals, this method could have broad applications.

Zoom In Zoom Out Reset image size

The self-diffusivity of bcc iron via the SIA dumbbell mechanism calculated using SEAKMC is compared directly with MD simulations in figure 10. To achieve desired statistics, the calculation of self-diffusion coefficient is based on multiple runs of 5000–10 000 KMC steps (successive jumps). For SIA and vacancy diffusion, one defect is introduced into a ∼2000 atom system; the defect concentration is 0.05%. Excellent agreement between SEAKMC and MD is observed at low and intermediate temperatures. Note that the SEAKMC does not require any a priori information on diffusion mechanisms or parameters; therefore, SEAKMC has the same predictive capability as MD. In addition, this agreement also confirms the accuracy of calculations in both migration energies and attempt frequencies. As temperature increases, SEAKMC starts to deviate from the MD results. This is because SEAKMC is based on harmonic transition state theory while MD follows the real dynamics. When the system is close to the minimum, such as low temperatures, the harmonic approximation adequately describes the potential well. As temperature increases, the anharmonicity of the atomic vibrations and the changes in potential energy landscape leads to a deviation from the harmonic approximation and the difference between SEAKMC and MD. The temperature at which such deviation starts depends on the interatomic potential and defect type. For instance, the SEAKMC results for SIA diffusion deviate from MD around 600 K while the agreement between SEAKMC and MD for vacancy diffusion holds up to ∼1000 K. Nevertheless, the difference between SEAKMC and MD is much smaller than one order of magnitude and is less than the standard error obtained in many MD simulations of diffusion at low temperatures (e.g. see the 700 K MD results in figure 11). The larger uncertainty in MD at low temperatures arises from the computational challenge with very long simulation times. Overall, this comparison clearly demonstrates that SEAKMC and MD are complementary to each other. At low temperatures, harmonic transition state theory is accurate and SEAKMC can be used to predict the dynamics with a much higher efficiency. At relatively high temperatures, the behavior starts to deviate from harmonic transition state theory and MD is more accurate. The behavior is similar for the defect diffusion coefficient (figure 10) and has been observed in previous studies [34]. As a large number of jumps were sampled in both SEAKMC and MD for SIA, the standard error is smaller than the symbol shown. The results of a similar study for vacancy diffusion in bcc iron are shown in figure 11. Compared with SIA, the difference between SEAKMC and MD reduced for all the temperatures calculated, indicating the vacancy has a more harmonic potential well than SIA. The effective energy barriers calculated from the Arrhenius equation for SIA and vacancy diffusion are 0.309 eV and 0.630 eV respectively, which is in excellent agreement with the previous calculations [24].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To further examine the accuracy of SEAKMC, the tracer correlation factor for the vacancy mechanism was calculated. The value obtained is 0.73, which is in excellent agreement with the theoretical value (0.727) [35]. The accuracy of the vacancy correlation factor demonstrates the fidelity of SEAKMC for atomistic processes. Therefore, we also calculated the tracer correlation factor for interstitial diffusion via the dumbbell mechanism, which is very time consuming to obtain using MD. SEAKMC predicts the tracer correlation factor for this case to be 0.44 if only the first nearest neighbor (FNN) jump mechanism is involved. The value would be higher if dumbbell rotation is included. The value depends on the relative probability for each process, which is determined by the saddle point energy. Interestingly, it is found that the FNN jump and rotation exhibit very different features since the FNN jump involves both defect diffusion and mass transport whereas rotation only involves mass transport. The prediction of the tracer correlation factor for the dumbbell mechanism demonstrates the unique predictive capabilities of SEAKMC.

In addition to diffusion, SEAKMC has also been applied to investigate defect interactions. For example, the recombination radius between an interstitial and vacancy is found to be ∼2.5 lattice constants using SEAKMC. Furthermore, the evolution of multiple defects was studied using an example observed in an MD displacement cascade simulation. The interactions between a nine-SIA cluster and an isolated vacancy and dumbbell SIA were studied using SEAKMC and shown to compare well with MD simulations [17]. The short-term comparison with MD confirmed the fidelity of the SEAKMC, while this method was able to extend the simulation time to a range that is well beyond the capability of MD.

In summary, SEAKMC provides a general framework for simulating defect evolution with atomistic details while requiring only the interatomic potential and no input of defect migration mechanisms, energies and corresponding attempt frequencies. SEAKMC consists of several components: AVs, saddle point searches, KMC and static relaxations. The method of determining the required AVs and the effect of AV size on saddle point energies were elucidated. A multi-step procedure involving the initial use of a small AV followed by successively larger AVs was proposed and has been observed to exhibit high computational efficiency. The observation that the dimer method has a higher probability of finding the saddle point with lowest energies is found to be associated with the symmetry of the defects. A general criterion to estimate the accuracy of the dynamics and simulation time was described and applied to representative scenarios. Self-diffusion due to SIA dumbbells and vacancies in bcc iron were used as applications to demonstrate the accuracy of SEAKMC. SEAKMC is in excellent agreement with MD at low to intermediate temperatures. At relatively high temperatures the loss of harmonic behavior leads to a deviation from MD. In addition, the tracer correlation factor for diffusion by the dumbbell mechanism is predicted and effects of FNN jump and rotation are discussed. Overall, the results demonstrate that SEAKMC can be used as a complementary tool to MD simulations. When harmonic transition state theory is applicable, SEAKMC has the same predictive capabilities as MD, but can reach much longer time scales.

Research sponsored by the US Department of Energy, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, 'Center for Defect Physics', an Energy Frontier Research Center.