Towards the reliable calculation of residence time for off-lattice kinetic Monte Carlo simulations

The simplest approaches to classifying kinetic events make use of individual macroscopic metrics such as the activation energy of an event and the relative energy of its adjacent minimum configuration [10], for example. However, these macroscopic values do not uniquely specify a given process, so researchers have turned to the use of more direct topological approaches.

Béland et al used a graphical isomorphism approach, specifically making using of a program called NAUTY [20, 21], to specify unique kinetic processes [22]. In this approach, a center-point for each event must first be identified. A sphere with radius, R, is then drawn around this center and the set of atoms that lie within this sphere then make up the nodes of a truncated graph. The edges of the graph are defined to connect any two nodes (atoms) that are closer than some edge cutoff distance, a, in the physical configuration. An example of how a graph is constructed in this fashion is shown in figure 2. A kinetic event, using this approach, is then completely specified by three graphs, one for each of the initial configuration, the transition state configuration, and the adjacent minimum configuration. NAUTY, or a similar tool, can be used to analyze these graphs and assign a set of hex keys corresponding to the canonical representation of each one. From this, the set of three hex keys that identify the graphs at the initial, transition, and final states uniquely specify a kinetic process, and anytime an event is identified, its uniqueness can be determined by comparing the hex keys associated with its configurations to those that have been found previously.

Zoom In Zoom Out Reset image size

This approach is powerful as it minimizes the amount of information that must be stored to fully specify a kinetic event and has a very straightforward comparison rule. However, it has a number of issues that make it ultimately unfavorable as a tool for uniquely specifying kinetic events. For example, the method is unable to account for the positional irregularities associated with the numerically converged configurations used in off-lattice methods. This results in either the false classification of two events as being distinct when they are, in fact, the same, or vice versa. In figure 2 an example is shown demonstrating how a small, numerically-based irregularity in an atom's position results in a graph with one additional edge and completely different graph hex keys. Additionally, the use of hex keys for transition state identification, though offering a simple means of defining kinetic processes, leaves the user no information with which to compare the similarity between two events that were found to be distinct, making it difficult to validate and optimize the method.

An alternate topological approach developed by Nandipati et al [23] that suffers from similar issues stores the transition state configuration as a binary number by overlaying a fine grid in the region central to the transition and assigning ones to the voxels of the grid that contain the center of mass of an atom and zeros to all other voxels. The grid is then unwrapped into a string of bits that make up a large binary number. This number can be conveniently stored and referenced. However, uncertainty of atomic position can, once again, confound this method, making it much too easy to store the same kinetic event multiple times due to minor perturbations in the atomic structure, complicating any attempt to determine the completeness of the energy landscape.

Due to concerns with the existing methods for uniquely specifying a kinetic process, we have pursued an alternate approach. Since our goal is to uniquely label atomic topologies without abstraction that would lose the physical features of the structure, we have chosen to work directly with atomic coordinates and trajectories. Furthermore, knowing the importance of uncertainty in atomic positions associated with numerical methods in general and with ART in particular, we hope to leverage that knowledge in event labeling. Specifically, we store atomic motion vectors associated with the transition. For each kinetic event found with ART, the x, y, and z displacements of each atom are stored between both the initial state and the saddle point, and the initial and final states. A schematic illustrating an event in which 4 atoms are displaced is shown in figure 3. In our method, the full vector description of each ${{t}_{i}}$ displacement is stored to uniquely specify the kinetic event.

Zoom In Zoom Out Reset image size

The comparison of two kinetic processes can be conducted by comparing the x, y, and z displacements of each atom in each of the events. If the difference between the displacements of any atoms in any direction is greater than a cutoff value $\epsilon$, which is chosen based on the known uncertainty associated with atomic position for the force tolerance being employed, then the events are taken to be unique. With a properly chosen value of $\epsilon$, atomic motions that are literally indistinguishable to within numerical uncertainty associated with the potentials and ART can be clearly identified; provided the uncertainty in atomic positions is deemed acceptable in the first place, the uncertainty in comparing and labelling events is guaranteed to be no worse.

At first glance, this approach might seem unnecessarily cumbersome compared with some of the memory-leaner alternatives previously discussed, but it has several crucial advantages. To begin with, it uses the most basic definition of a kinetic event, i.e. the movements of atoms in the system, to define and distinguish kinetic events. Secondly, because it makes use of the known numerical uncertainty of atomic positions in distinguishing between events, it allows the user to choose convergence criteria that are logically consistent with an application without sacrificing the ability to compare events. A comparison of the present method with previous approaches is presented later in the Implementation section of this paper.

If the information about each kinetic process is stored in a database, we can easily minimize the search space by performing comparisons with only the subset of events that have, e.g. activation energies within some tolerance of that found in the current search, or some other rapid screening criterion. This has the added advantage of allowing us to consider only the displacements between the initial and final configurations of a process, which halves the number of required comparisons by excluding the saddle point, which is associated with greater uncertainty in atomic position. Though it is possible for multiple activated states to lead to the same final state, we have not found these activated states to have similar enough activation energies so as to be confused with one another, at least in the system of our analysis. It is possible that some systems may be able to access symmetrically related transition states, in which case, both the saddle point and final configurations would need to be compared to identify a unique event.

The simplest approach to terminating an ART search is to qualitatively deem the search to be complete after some arbitrary number of search attempts or successful searches [17, 18]. For example one could deem their search complete when n unique saddle points have been identified [17, 24] but n remains somewhat arbitrary. More sophisticated stopping criteria have been introduced; for example Fan et al in [24] terminate a search if a newly found event has a kinetic probability that is less than 0.1% of the current residence time estimate. However, this approach still relies on arbitrary factors such as the order in which events are randomly identified. An improvement over arbitrary termination is to choose some property of the search or the energy landscape around which a statistical test can be built to determine the probability that the chosen property is sufficiently close to the true value. For example, in [10] we used the frequency of the most commonly found kinetic event as a metric for terminating the search. Specifically, we binned all the kinetic events found with random ART searches and tracked how the average frequency across the bins of the most commonly found saddle point, ${{f}_{\max}}$, changed with additional successful searches. We then used the sample mean of ${{f}_{\max}}$ and its standard deviation, s, in a Student's t-test to calculate the minimum sample size, ${{n}_{\min}}$, that would be required to achieve the desired confidence that our calculated value was sufficiently close to the true value, and once our sample size exceeded ${{n}_{\min}}$, we deemed our search complete.

The above proposal was presented in the context of a static energy landscape search, i.e. without connection to a KMC context or, indeed, any structural evolution. Since system kinetics are of prime interest to a KMC application, we presently suggest that a further improved approach would be to define a metric based on a kinetic property of the system, rather than a numerical one, as has been previously suggested [24]. Specifically, since an accurate calculation of the residence time, ${{\tau}_{\text{R}}}$, is our goal, we propose to use the convergence of ${{\tau}_{\text{R}}}$ to assess the completeness of an energy landscape search. The metric used to establish convergence in our present approach is discussed below.

Whether one uses an arbitrary or quantitative approach to search termination, the random nature of the perturbations typically used in the search method are a common limitation. In either case, we cannot be sure that we have effectively searched the energy landscape, and, in either case, we must perform a large number of searches that, by definition, contain no new information in order to sufficiently characterize the energy landscape; this all results in an inefficient use of computational resources.

As an alternative to the stochastic search method, we employ a deterministic search method by applying a set of systematic perturbations (each of which separately initiates the ART algorithm) to one atom in the system at a time. The perturbations make up the points of a 3D mesh, with points that are uniformly spaced on concentric spheres with radii corresponding to the magnitude of the perturbation, ${{d}_{i}}$, and centered on the atom being investigated. The product of the number of directions on a perturbation sphere, $\rho$, and the number of shells in the radial direction of the mesh, $d$, is the total number of perturbations in the search of each atom. In figure 4, perturbation meshes with several values for $\rho$ and a representative set of concentric shells (with $d=4$) that are separated in the radial direction by a distance δ, are shown to illustrate the nature of the perturbation space used in the deterministic search.

Zoom In Zoom Out Reset image size

When $\rho$ and d are made sufficiently large and $\delta ~$ is made sufficiently small, the systematic search can explore all the relevant kinetic pathways in a system. However, determining the critical values for $\rho$, d, and $\delta$ requires some optimization. By calculating the residence time for a system as a function of these parameters, a critical value—beyond which residence time is independent of $\rho$, d, and $\delta$ respectively—can be determined for each of them. One caveat important to the validity of this approach is that meshes of increasing $\rho$ or d must not inherit all the same perturbation directions as the coarser meshes. This requirement removes the bias on $\rho$ and d that is introduced if particular directions or radii in the coarse mesh are very advantageously aligned for searching certain atomic environments.

It should be noted that future approaches may achieve efficiency improvements in this search process by exploiting the symmetry in a LAE to minimize the set of perturbations required to search its energy landscape. Additionally, as with other methods that rely on the perturbation of individual atoms to initiate the search for kinetic events, there is the necessary concern that concerted events involving large numbers of atoms will not be identified by the search method. In our prior work, we have found that eigenvector-following methods such as ART allow for collective mechanisms to be identified due to the inherently non-local nature of the eigenmodes of a system [10]. This does not guarantee that all collective mechanisms will be found during a local perturbation search process, as described in this paper. As such, future studies may pursue a more exhaustive investigation into this topic.

We consider values for $\epsilon \,$ in the range $0.010-0.050\,{\mathring{\rm A}}$ based on our knowledge that atomic position has an uncertainty on the order of $0.01\,{\mathring{\rm A}}$, but can be as much as $0.025\,{\mathring{\rm A}}$ in minimum energy configurations. To assess how the residence time is affected by $\epsilon$, we have calculated the residence time, according to equation (1), using our systematic search method for various choices of $\epsilon$ in this target range. Since the critical value of $\rho$ is initially unknown, we perform the systematic search for a range of $\rho$ values (with d  =  4, $\delta ~$  =  0.5 Å, and ${{d}_{1}}=0.75\,{\mathring{\rm A}}$ held constant), the results of which are shown in figure 6.

Zoom In Zoom Out Reset image size

From figure 6, we note that the calculated residence time seems to be relatively converged for about 500 perturbation directions; a more quantitative assessment of this is discussed below. As expected, larger values of $\epsilon$ result in larger residence times due to cases where distinct events are incorrectly identified as being the same, while the smallest values of $\epsilon$ result in the opposite effect, where the same kinetic event is incorrectly counted multiple times. The optimal value of $\epsilon$ is different for each event. Some kinetic events are somewhat sloppy, with the highest uncertainty in their atomic positions ($\pm 0.025\,{\mathring{\rm A}}$), for which a large $\epsilon$ is preferable, while many other events are much cleaner and a smaller $\epsilon$ could be used. However, for meshes with greater than 250 perturbation directions, the difference between the residence time calculated with the different choices of $\epsilon$ is, at most $7.9\centerdot {{10}^{-11}}\,~\text{s}$. Thus, we take half this value, $4.0\centerdot {{10}^{-11}}\,\text{s}$, to be the uncertainty associated with the residence time we obtain with this method (shown as error bars in figure 6) and choose an intermediate $\epsilon =0.025\,{\mathring{\rm A}}$ in future calculations.

At this point, it is important to recall that none of the exact points on the sparser meshes are inherently included in the denser meshes. Thus, it is possible to find certain kinetic events with a sparser mesh that are then missed with the denser mesh. This prevents the determination of the minimal mesh size from being biased by a 'lucky' sparse mesh population that is inherited by all denser meshes and accounts for the non-monotonic behavior seen in figure 6.

Using $\epsilon =0.025\,{\mathring{\rm A}}$, we can now optimize the mesh parameters required for a complete search of the energy landscape for our target system. From figure 6 it appears that the critical value of ρ may be around 150; however, it is possible that with different choices of d and δ a more minimal perturbation sphere may be used. We proceed by optimizing the radial parameters of the mesh. After finding that no search converged to a saddle point when the initial perturbation was less than 0.75 Å in our implementation of ART, we define the space of possible radial perturbation distances as ranging from 0.75–2.25 Å. We consider shell spacing values of $\delta$ ranging between 0.1–1.0 ${\mathring{\rm A}}$. Setting $\rho =328$, we investigate residence time as a function of d and $\delta$, as shown in figure 7.

Zoom In Zoom Out Reset image size

From figure 7, we note that the residence time is independent of d so long as $d~>~4$ and $\delta =0.15~\,{\mathring{\rm A}}$. As such, we have found that including perturbations $>1.20~\,{\mathring{\rm A}}$ did not affect the residence time calculation. Decreasing $\delta$ to 0.10 Å did not affect the residence time either. However, we note that the residence time did show sensitivity to the choice of the lowest ${{d}_{i}}$ value, ${{d}_{1}}$. Holding d  =  5 and $\delta =0.15\,{\mathring{\rm A}}$, we increased ${{d}_{1}}$ from 0.75 ${\mathring{\rm A}}$ to 0.85 ${\mathring{\rm A}}$ and observed an increase in residence time.

Having a sense of the optimal set of perturbation shells with ρ held constant, we next need to determine the codependence of our mesh parameters. In figure 8, we show the residence time as a function of d for $\rho ~=~156$ and $\rho ~=~328$ with $\delta ~=~0.15$ ${\mathring{\rm A}}$ and ${{d}_{1}}=0.75\,{\mathring{\rm A}}$ and error bars centered on data with $\rho ~=~328$. In figure 8, it can be seen that the residence time versus d behavior is independent of the mesh density ρ, with both datasets converging to d  =  4 as the minimal number of radial shells in the perturbation mesh.

Zoom In Zoom Out Reset image size

Having established that the mesh parameters are independent, we can now revisit figure 6 and determine the minimal mesh density, $\rho$, for this system. Noting that choosing a constant value for $\epsilon$ lead to an uncertainty in the residence time of $4.0\centerdot {{10}^{-11}}\,\,\text{s}$, we can calculate the standard deviation of the residence time for the data shown in figure 6 with $\epsilon =0.025\,{\mathring{\rm A}}$. By successively neglecting the data from the coarsest meshes, we obtain a series of standard deviation values, ${{s}_{\overline{{{\tau}_{\text{R}}}},i}}$, where ${{\overline{{{\tau}_{\text{R}}}}}_{,i}}$ is defined as the average residence time calculated using the set of meshes with $\rho \geqslant ~{{\rho}_{i}}$ and ${{\rho}_{P}}$ being the finest mesh:

$$\begin{eqnarray}&&{{\overline{{{\tau}_{\text{R}}}}}_{,i}}=\frac{1}{P-i}\underset{j=i}{\overset{P}{\sum}}\,{{\tau}_{\text{R},{{\rho}_{j}}}}\end{eqnarray} \tag{ 4 }$$

We use the requirement that a minimally sufficient mesh must satisfy ${{s}_{\overline{{{\tau}_{\text{R}}}},i}}<4.0\centerdot {{10}^{-11}}\,\text{s}$. In figure 9 we plot ${{s}_{\overline{{{\tau}_{\text{R}}}},i}}$ as a function of $\rho$ (along with ${{\tau}_{\text{R}}}$ for illustrative purposes). The uncertainty associated with ${{\tau}_{\text{R}}}$ is shown as a horizontal black line; the point at which ${{s}_{\overline{{{\tau}_{\text{R}}}},i}}$ crosses below this line, at $\rho =46$, is marked with a vertical dotted line, and we can take $\rho =46$ to be minimally sufficient for the present system of study.

Zoom In Zoom Out Reset image size

In figure 10, we show several examples in which we have used NAUTY to identify unique kinetic events and LAEs for the energy landscape search of the  ∑5(2 1 0) grain boundary for several example combinations of the R and a parameters (recall that R is the radius of the truncated graph and a is the cutoff distance for an edge). For the purposes of this comparison, we have only used meshes with ρ in the most relevant range: 46–562. To allow for a direct comparison with the results of figure 6, we have held d  =  4, $\delta =~0.5$ ${\mathring{\rm A}}$, and ${{d}_{1}}=0.75\,{\mathring{\rm A}}$ constant.

Zoom In Zoom Out Reset image size

From this figure we see, as was observed with our method above, that the residence time appears to converge for perturbation meshes with greater than about 500 directions. The sensitivity of the residence time to the choices of R and a are notable, as the spread of calculated residence times in the converged region, ~$1\centerdot {{10}^{-10}}$ s, is about the same as the residence time itself, and there is no obvious way of determining which combination of R and a lead to the true residence time of this system.

In figure 11, we plot the same data as shown in figure 10 and compare it with the data from figure 6 with $\epsilon =0.025\,\,{\mathring{\rm A}}$. From figure 11, we note that residence time calculated using the present method lies within the same range as that calculated using NAUTY, but with a tighter uncertainty (under half). In all cases, the meshes with greater than 500 directions show a converged residence time behavior.

Zoom In Zoom Out Reset image size

In addition to having a tighter uncertainty, our approach has the advantage of not depending on the arbitrary choice of parameters R and a. Because it relies instead on the physical definition of a kinetic event (the displacement of atoms) in distinguishing between events, we believe it provides a reasonable means of determining the residence time in an off-lattice atomic configuration.

To determine the minimum number of closest neighbor atoms, M, required to specify an LAE for our system, we have calculated the residence time of our example system using values of M from 12 to 64. Only the atoms in or near to the grain boundary were of interest in the kinetic analysis, thus atoms more than one nearest neighbor spacing from the GB were not searched for kinetic events. In this highly ordered GB, the number of unique LAEs in the GB region ranged from 8 to 22 for the investigated values of M. The results are shown in figure 12. The error bars denote the uncertainty in the residence time calculation. Due to the minimal fluctuation observed in the residence time across the full range of LAE sizes considered, we conclude that 12 closest neighbors is sufficient for the definition of the LAE in this case.

Zoom In Zoom Out Reset image size

In figure 13 we show the perturbation directions that resulted in successfully finding a saddle point for three different LAEs for perturbation meshes with ρ  =  128 and ρ  =  562. Perturbations that led to a successful ART search are plotted on a unit sphere, or 'transition globe', and colored according to the activation energy of the corresponding kinetic event. For clarity in visualization, perturbations of all magnitudes are collapsed onto a single globe.

Zoom In Zoom Out Reset image size

From figure 13, we see that in the cases where a denser mesh was used, there are 'continents' of perturbation directions that result in what appear to be the same kinetic event. In some cases, these continents are quite large, indicating that a broad range of perturbations will lead to convergence to the same saddle point. In other cases, only a narrow set of directions lead to a given event. It is these latter kinetic events, if sufficiently low in activation energy, that dictate the required mesh density of a system and which highlight the risks of declaring any search 'complete': there is always some chance of missing the critical perturbation that might find a low energy (highly kinetically significant) transition.

To investigate the potential relationship between the activation energy of a kinetic event and the ease with which it can be converged upon, in figure 14 we plot the solid angle, in steradians, that a 'continent' on a transition globe occupies as a function of its activation energy. In general, we see that across the whole range of activation energies there are many events that will only be observed if triggered by a relatively localized perturbation. Conversely, there is a subset of transitions for which a broad range of perturbations will converge to what appears to be the same saddle point, as can be observed in figure 13. There is a possible correlation between higher energy saddle points with larger continent areas. One explanation of this trend might be that higher energy saddle points are less likely to fall back down into the initial minimum state during convergence than lower energy activated states, so a wider set of initial perturbations can converge to the same saddle point in these cases.

Zoom In Zoom Out Reset image size