Long time scale simulation of a grain boundary in copper

A general, twisted and tilted, grain boundary in copper has been simulated using the adaptive kinetic Monte Carlo method to study the atomistic structure of the non-crystalline region and the mechanism of annealing events that occur at low temperature. The simulated time interval spanned 67 μs at 135 K. Similar final configurations were obtained starting from different initial structures: (i) by bringing the two grains into contact without any intermediate layer, and (ii) by inserting an amorphous region between the grains. The results obtained were analyzed with a radial distribution function and a common neighbor analysis. Annealing events leading to lowering of the energy typically involved concerted displacement of several atoms—even as many as 10 atoms displaced by more than half an Ångström. Increased local icosahedral ordering is observed in the boundary layer, but local HCP coordination was also observed. In the final low-energy configurations, the thickness of the region separating the crystalline grains corresponds to just one atomic layer, in good agreement with reported experimental observations. The simulated system consists of 1307 atoms and atomic interactions were described using effective medium theory.


Introduction
The structure and properties of grain boundaries (GBs) in metals have been extensively studied as they can have a key influence on the properties of the material. A large number of papers using atomic-scale computer simulations to study GB properties-or properties controlled by GBshave therefore appeared in the literature in recent years. For example, molecular dynamics (MD) simulations have shown how the presence of GBs in nanocrystalline metals leads to a maximum in hardness at a specific grain size [1]- [3] and how the GBs influence the nucleation of dislocations in nanocrystalline metals [4]- [6]. Trijunctions where GBs meet have also been simulated and the energetics and atomic structure have been studied [7].
The nanometer grain size makes it possible to simulate a complete grain structure of a nanocrystalline metal with MD and the role of high-temperature relaxation in this domain has been analyzed [8]. An alternative approach is to focus on a single flat GB in order to study its properties. MD studies of such systems have been used to study self-diffusion in GBs [9,10], GB migration [10]- [12] and recrystallization [13].
Common to all these applications of MD to GBs is the timescale problem: MD can only describe processes occurring at timescales of at most a few nanoseconds. To alleviate this problem, many MD simulations operate at high temperature or with high driving forces for the processes being studied. For the timescale, a number of accelerated dynamics schemes have been proposed, including hyperdynamics [14], the parallel replica method [15], temperature accelerated dynamics [16] and adaptive kinetic Monte Carlo (AKMC) [17]. All of these methods operate under the assumption of infrequent events, i.e. that the system can be described as 3 undergoing a series of discrete transitions between more or less well-defined states (see e.g. the review by Voter et al [18]).
In this paper, we demonstrate that infrequent-event dynamics can be applied to a materials science problem, namely the annealing of a GB at low temperatures, using the AKMC algorithm. The simulations use two very different starting configurations to reach similar final structures for the GB, where the GB is very narrow.
The atomic-scale structure of GBs in nanocrystalline metals has been debated since the suggestion that GBs in nanocrystalline metals have a 'frozen gas' structure with a low density [19]. More recent experimental studies have contradicted this, indicating that the thickness of the GB region is limited to 0.4-0.5 nm, leaving only a few atomic layers without any crystalline structure as in coarser grained materials [20]- [22]. MD simulations of nanocrystalline metals typically show similarly narrow GB structures, although care should be taken when interpreting this as evidence for narrow GBs in nanocrystalline metals, since the initial configurations for these systems by construction have narrow GBs [23,24]. The few MD simulations of a solidifying nanocrystalline metal also yield narrow GBs [25]. In this paper, we confirm this by simulations over a much longer timescale (up to 67 µs), and show that annealing events occur during the entire period, leading to lower energy GB structure than what could have been found with conventional MD.

AKMC
The AKMC method is an extension of kinetic Monte Carlo (KMC). Here, we apply AKMC to obtain atomistic dynamics over long timescales. The adaptive label refers to the creation of a table of events (TOE) on-the-fly during the AKMC dynamics without a priori knowledge. To find reaction mechanisms, we use a minimum-mode following method [26,27] to locate saddle points (SPs) on the potential energy surface (PES). When located, the escape rate is calculated and the mechanism added to the (adaptive) TOE [17]. Only those mechanisms determined by the SP searches are included in the simulation. For many systems, however, one knows certain characteristics of the low-barrier processes. For example, ad-atoms diffusing on a surface can be characterized as active because they have low coordination numbers. This kind of knowledge can be used to improve sampling in the active regions of configuration space with a limited number of SP searches.
When searching for an SP, the force acting on the system is modified according to the eigenmode of the Hessian with the lowest eigenvalue. If this lowest eigenvalue is positive, the force component orthogonal to the eigenmode is zeroed out and the parallel component is reversed. In other words, the eigenmode is followed up the PES. If a negative eigenvalue is obtained, the forces are modified by reversing the component parallel to the displacement in the eigenmode, allowing the system to relax along all perpendicular directions. This modified (effective) force for λ low,i 0, is followed using a conjugated gradients method to reach an SP. The lowest eigenvalue mode is found with the dimer method [28]. Sampling of mechanisms from the reactant configuration is done with a stochastic displacement from the minimum to initiate each SP search.

4
The rate constant of a mechanism is calculated by assuming that the PES is harmonic at the SP and minima. Under this assumption, harmonic transition state theory (HTST) applies and the rate for passing through a specific channel is calculated as where E S is the energy of the SP, E R is the energy of the reactant configuration, ν S,i and ν R,i are the frequencies of the vibrational modes at the SP and at the reactant configuration, respectively and D is the dimensionality of the Hessian. When an escape channel is located, its rate constant and product state are determined and tabulated in the TOE. When the table is considered complete, a process is chosen by the KMC algorithm and the corresponding time step calculated as where µ is a random number in the interval (0, 1) and j runs over all distinct mechanisms.

Distributed SP searches
The task of creating the TOE is by far the most computationally demanding part of an AKMC calculation. These searches are independent, however, and can be done in parallel on a cluster or network of distributed computers. The current simulation was carried out using EON 5 , which is a distributed system based on a master-slave (server-client) paradigm. During an AKMC iteration, the master continuously sends displaced configurations to nodes in a public farm of slaves. The slaves then perform SP searches and report the results back to the master where they are tabulated. When the master has a sufficiently converged TOE, a mechanism is selected with a probability proportional to its rate and the time step is determined by equation (4). The system is then advanced to the final state of the chosen process where the algorithm is repeated [29].

AKMC parameters
The focus of this work is the atomistic dynamical behavior at the GB. To sample the relevant parts of configuration space, the epicenter for the displacement used to initiate an SP search is restricted to the subset of atoms that belong to the GB. The GB atoms are defined as those that do not have face-centered cubic (FCC) coordination. A displacement vector is created including all atoms that are within a sphere of radius 6 Å from the epicenter, regardless of their coordination. The components for the displacement vector are drawn from a Gaussian distribution with a standard deviation of 0.3 Å. The initial configuration for the AKMC annealing simulation consists mainly of crystalline ordered FCC atoms (see figure 6(A3) or table 3). In the initial quenched configuration, the total number of GB atoms is about 200. In each AKMC step, 200 SP searches are assumed to be sufficient to find the significant reaction mechanisms. To validate this assumption, two Table 1. Statistical distribution of 1000 SP searches. Two configurations have been analyzed (*1 and *2). Mechanisms labeled A-C are the three lowest barrier events, whereas D is the one that is most frequently determined. The column labeled 'Moved' accounts for atoms that are displaced by more than 0.5 Å in the SP configurations. The full data series for system 1 is shown in figure 1.  The mechanism labeled A involves two atoms ( 0.5 Å) and accounts for 3.1%. B involves 4 atoms and accounts for 1.7%, C involves four atoms and accounts for 0.9%. The mechanism labeled D is marked due to its abundance; it involves four atoms and accounts for 7.3%. As the lowest energy mechanism is determined once for every 33 searches, the criterion of 200 converged SPs should be more than sufficient to ensure a determination of the most significant mechanism.
configurations from the simulation were selected for further analysis. For each of the two configurations, 1000 SP searches were performed. The distribution of the SPs is shown in table 1 and figure 1. The lowest barrier process from the two states is found on average every 33 and Table 2. EMT parameters for copper, in units derived from eV and Å. The parameters are the default values as defined in [31].   19 searches, respectively, so that our requirement of 200 converged SPs appears to be sufficient in most cases. On average, all of the low-barrier mechanisms tabulated would be accounted for within 200 searches.

Interaction potential
The copper interactions are described with an effective medium theory (EMT) interatomic potential [30,31] using the parameter set in table 2. The EMT potential has been used for more than a decade and this parameter set for copper has successfully been used to describe the behavior of clusters, surfaces and solid structures [1,32,33].

Structural setup and preparation
To synthesize a single twisted and tilted general GB, two misaligned single-crystalline domains are brought into contact. Periodic boundary conditions apply in theX -andŶ -directions that span the GB plane; theẐ -direction is normal to the boundary. The two outermost layers in thê Z -direction are kept frozen at bulk FCC lattice sites to minimize surface effects. The structure consists of 1309 atoms, which corresponds to 22 layers of 60 atoms each in an ordered FCC crystal.
For the simulation, two different initial structures were generated. A thin GB was modeled by placing the atoms at FCC lattice sites in both the upper and lower domains. The interface itself is between the two perfect surfaces. The other wide GB has six interface layers that are initially amorphous. Schematics of the initial structures are shown in figure 2. It should be noted that the two initial structures with thin and wide boundaries are thinner and wider (respectively) than the structures modeled with AKMC due to an initial annealing step.

Thin GB
The thin GB consists of two single-crystalline FCC grains that meet in a general GB. In principle, any GB modeled with periodic boundary conditions is a coincidence site lattice GB, but in this case the ratio of coinciding lattice sites is very low ( = 8580) and it is more reasonable to regard it as a general GB. Both crystalline domains have the sameX -and Y -periodicity in the supercell and periodic boundary conditions apply in these two directions. The GB is acquired by faceting the two grains to obtain surfaces with the supercell'sẐ -direction as their normals. These two parallel surfaces are then brought into contact. The initial distance between the closest atoms in the lowermost and uppermost domain is 1 Å. From the insets in figure 12, it can be seen that the interfacing surfaces are different. The interface above the boundary could be considered an FCC (001) surface and that below an FCC (111) surface. The thin GB initially contains 161 GB atoms and is shown in figure 3(A1).  ) and (B2) have similar distributions between the two initial configurations, which indicates that a similar annealed GB structure has been obtained from both our thin and wide initial structures.

Wide GB
The wide GB is constructed from a thin boundary sample. A slab of approximately six layers is sliced out from the center region of the sample, which is melted with Langevin dynamics at 1500 K for 10 ps. The liquid slab is then quenched and the resulting amorphous structure is reinserted into the interface. The wide GB initially contains 519 GB atoms and is shown in figure 3(B1). This subset of GB atoms is used to obtain the radial distribution functions (RDFs) in our later analysis.

Structure preparation
The initial artificial structures were annealed with MD using a Langevin thermostat at 300 K for 250 ps. This treatment increases the number of GB atoms for the thin boundary to 267 and decreases it to 261 for the wide. The structures are further annealed for 13 ns at a lower temperature of 135 K. The Z-coordinates were then rescaled to reduce any strain induced by the frozen outermost atoms. The configurations before and after annealing are depicted in figure 3 and their RDFs are shown in figure 4. Table 3 summarizes the evolution of GB atoms during the preparation processes.
The thin GB configuration was selected to undergo AKMC long timescale annealing, although we expect similar results for the initially wide boundary because their initial structures are so similar.

Results
The output from an AKMC simulation represents the distinct configurations visited during the simulation and a list describing in which order and for how long these configurations were visited. This enables one to evaluate both state properties and detailed atomistic analysis of specific transition mechanisms. To analyze state properties, we focused on interface energy and two general order properties: the RDF and common neighbor analysis (CNA). These analysis tools were used to characterize the mechanism of two events that are believed to be of particular importance. The interface energy is defined as where E conf is the structural energy of the configuration, E surf is the energy of the two surfaces facing vacuum, E FCC is the energy of a perfect FCC crystal with the same number of atoms as the simulated structure and A is the footprint of the configuration. As a reference, the energy of a relaxed copper (001) FCC surface facing vacuum is 69.2 meV Å −2 .

Interface energy versus time analysis
We have used two ways to visualize the interface energy, as a function of AKMC iterations (figure 5 inset) and as a function of accumulated time ( figure 5). In the energy versus time graph, a long-lived configuration shows up as a plateau; rapid transitions appear as spikes. In the time versus energy plot (figure 6), two significant energy lowering events are apparent. The first occurs at 55 µs and causes an energy drop of 0.24 eV, corresponding to an interfacial energy lowering of 0.6 meV Å −2 . The second, which happens after 66 µs, causes the energy to lower by 0.08 eV and the interface energy by 0.2 meV Å −2 . The atomic rearrangements during these two events are subject to a detailed analysis in sections 4.4.1 and 4.4.2. A general effect of both events is to increase the frequency of transitions after they have taken place. It is remarkable that the system starts with transformations through a series of high-energy configurations before the system manages to reach a region of low-energy configurations, figure 6; a similar behavior has also been observed in earlier test simulations. The highest configuration encountered is 0.48 eV above the initial structure. Table 3 states the changes in the number of GB atoms from the initial structures to the fully annealed configurations. The annealed configuration at 300 K contains approximately 265 GB atoms, and at 135 K an optimum structure contains approximately 240 (A2 and B2). When the thermal distortion is removed and the structural energy solely defines the favorable structure, the number of GB atoms decreases further to 210 (A3 and B3). The dependence of the interface energy on the number of GB atoms for the sampled configurations during the AKMC annealing is depicted in figure 7. It appears that configurations containing between 198 and 212 GB atoms have the lowest interface energy. A monolayer consists of 60 atoms, which means that configurations where the interface has three layers without crystalline order are energetically favorable. Two other striking features are that all configurations containing a single HCP-coordinated atom have interface energies below 82.5 meV Å −2 . Furthermore, additional test simulations also evolve into configurations containing a single HCP atom. Thus, a locally HCP-coordinated atom appears to be a characteristic of the GB structure for this system size.

Structural analysis
The structural change during the AKMC simulation has been characterized with RDFs in figure 8. The main difference between the distributions before and after AKMC annealing is the sharper peaks, corresponding to an increase in crystallinity.
CNA [34] is conducted for all configurations to classify the local crystalline structure according to the distribution of common bonds shared by nearest-neighbor atoms. Important numbers are CNA421, CNA422 and CNA555 which indicate FCC, HCP and icosahedral structures, respectively. A cutoff radius of 3.3 Å is used to define neighboring pairs. The evolution of selected CNA values is shown in figure 9. The data are smoothed using a Gaussian weighted average of the 50 iterations before and after the actual iteration. The selected curves show the most abundant CNA values except for CNA421, since it corresponds to the bulk FCC coordination. Oscillations in the plateaus of the CNA curves are in phase with fluctuations in the interface energy, with the exception of CNA422 and CNA433. Sharp positive spikes are also in phase with the spikes in energy. Again CNA422 is an exception to this rule; the outof-phase behavior between energy and CNA422 is a clear indicator that it is a favorable local structure. At the energy lowering event occurring after 55 µs, one GB atom obtains a local HCP coordination, raising the CNA422 and CNA555 curves. After this event, the amplitude of fluctuations for all curves decreases, indicating that the configurations sampled are more alike.
The CNA data for a single configuration were analyzed to find the spatial distribution of a specific CNA value. The analysis was carried out by slicing a configuration into 31 equally sized slabs along theẐ -axis to define bins that the CNA values are sorted into. CNA555, indicating icosahedral structure, has a clear trend. The data for the initial and final configurations of the AKMC simulation are shown in figure 10. During the AKMC simulation, there is a build up of CNA555 pairs. The position of the maximal values and the tails of the two curves coincide, which means that changes have taken place inside the GB region. The maximal value is in the middle of the structure, where an increase of approximately 50% takes place during annealing. The tails above and below the boundary show substantial differences. The tail below increases  The result of the AKMC annealing is that the spikes and valleys become more distinct, indicating a transformation toward an ordered FCC structure.
abruptly, which indicates a flat interface between the single crystal and the GB region whereas the tail above is wide, which indicates that the roughness on this side is larger. It is in the tail above where the largest build up of CNA555-coordinated atoms occurs.  MD-annealing (A3) Figure 10. Layered CNA for CNA555 pairs tracing out the evolution of local icosahedral structure. A number of such pairs are partitioned into 31 layers along theẐ -axis to form a histogram. The area under these curves is 30 for the MDannealed structure and 60 for the AKMC-annealed structure. It is remarkable that the curves rise most rapidly when the boundary is approached from below. The plot shows that changes toward CNA555 coordination occur primarily in the upper part of the interface.

Event analysis
The two energy-lowering events are analyzed here in more detail. The first, at 55 µs, forms an HCP-coordinated atom and the second, at 66 µs flattens the GB.

Event at 55 µs.
When compared with the initial configuration for the AKMC simulation, the first event lowers the energy by 0.24 eV and the interfacial energy from 83.2 to 82.6 meV Å −2 . The curves for interface energy and CNA422 are shown in figure 11, where it should be noted that the data are depicted as a function of AKMC transitions and comparison with figure 6 should be done with care. The most striking features are the sudden increases on the CNA curve. At step 100, both the energy and the CNA curves increase. An out-of-phase change occurs at 500 transitions where the energy decreases and the CNA422 value increases. At this point, a single GB atom obtains HCP coordination. Since copper is an FCC metal, the HCP-coordinated atom can be considered a stacking fault. The transition creating the HCPcoordinated atom is a collective rearrangement where five atoms are displaced by more than 0.5 Å. The inset in figure 11 shows this collective rearrangement. Surprisingly, the most active A remarkable feature of the process is that an atom obtains a local HCP coordination. The transition where the HCP core is created involves five atoms being displaced by more than 0.5 Å. The core atom stays HCP coordinated throughout the rest of the simulation. The inset shows the creation mechanism where GB atoms are removed; exceptions are the five atoms that moved more than 0.5 Å, which are colored blue, and the GB atom that became HCP coordinated, which is green.
atoms are not the nearest neighbors to the formed HCP-coordinated atom. The energy is lowered due to structural relaxation around the newly formed HCP core. After the first event, the amplitude of fluctuations of all analyzed structural properties, energy included, decreases significantly, and the structures are more alike. The HCP core appears to be a very stable structure since it and its nearest neighbors maintain the local HCP coordination throughout the remainder of the simulation, if some minor short intervals are ignored.

4.4.2.
Event at 66 µs. The second event lowers the energy by 0.08 eV and the interface energy from 81.9 to 81.7 meV Å −2 . Figure 12 shows the energy as a function of AKMC transition. On the energy curve, there is only one significant step. This step accounts for 20 AKMC transitions and the corresponding structural change is a flattening of the GB. The changes in coordination caused by the flattening events are tabulated in table 4, where it should be noted that 10 GB atoms obtain FCC coordination.
The insets in figure 12 show that the interfaces between the GB and the two perfect singlecrystalline domains are different. The difference is significant since the activity on the (001)-like  Figure 12. Energy as a function of AKMC transitions during the second significant annealing event. The energy curve has been smoothed by a weighted Gaussian average of the 20 transitions before and after the actual transition. The total process caused a lowering of structural energy by 0.08 eV, changing the interface energy from 81.9 to 81.7 meV Å −2 , and took place over 20 transitions. The insets show theXŶ -planes of the GB as seen from above and below. Only GB atoms are shown, gray atoms were GB atoms both before and after the process, black atoms became GB atoms and white atoms obtained FCC coordination. The grids are imposed to help show that the interface below is a (111)-like surface, whereas above it is a (001)-like surface. In general, the reshaping is a flattening of the GB. surface above is more pronounced than on the (111)-like surface below. This is related to the observation (see figure 10) that the interface at the (001) high-energy surface is rough and that of the (111) low-energy surface is sharp.

Algorithm performance
The AKMC algorithm [17] has previously been used to describe surface growth and diffusion [35] where high prefactors were found to be important. In contrast, in this disordered bulk environment, 90% of all processes have prefactors between 10 13 and 10 14 Hz, and there is no apparent correlation between the individual processes' energy barrier and the prefactor (no compensation effect). We have also investigated the abundance of concerted displacements. In our analysis, an atom is considered displaced if it moves by more than 0.5 Å during a transition. A histogram of the number of displaced atoms is shown in figure 14. Transitions involving four or six atoms are nearly as frequent as transitions involving one or two atoms. This is a clear indication that KMC simulations that rely on a predefined TOE should be used with extreme care, as the a priori knowledge of these concerted mechanisms is unlikely.
In our simulation, the height of the energy barriers is often less than 0.1 eV. These low barriers have an unfortunate impact on the achievable timescales because each time step is dominated by the lowest barrier process: the highest rate from equation (4). Furthermore, in many cases the low-barrier events connect a limited subset of configurations. When the system enters one of these subsets, it is only able to leave again by overcoming a relatively high barrier, so that the system rattles between configurations that are revisited over and over again. In general, this problem is handled by storing the TOE and reusing it when the configuration is revisited. However, the average barrier height decreases during the simulation, and as a consequence the ratio of rates within a subset and leaving the subset increases. At a certain point, this causes the rattling problem to become so severe that the limiting factor for time evolution is how fast the KMC transitions can be performed. The simulation was terminated at this point. Figure 13 shows how time evolves along the simulation; in the first 25% (50 000) of the AKMC transitions, 75% (67 µs) of the simulated time has elapsed. Histogram showing the distribution of atoms being displaced by more than 0.5 Å. The analysis extends over the first 400 000 transitions. When the four most abundant transition types are accounted for (1, 2, 4 and 6 atoms), only 4% of the transitions remains. A transition involving either four or six atoms occurs almost as frequently as a transition involving only a single atom.

Discussion and summary
A GB in copper has been subject to AKMC annealing at 135 K for 67 µs. In the final configuration, the GB consists of approximately three atomic layers of atoms. The upper and lower interfaces to the single-crystalline domains have a high degree of order and it is only one of the three layers that is without any FCC coordination. The remaining two layers can be classified as the surfaces of the crystalline domains. The final GB thickness of 0.5 nm is in excellent agreement with experimental data [21,22]. When compared with earlier simulations, our results are similar to what Keblinski et al reported in that the ordering of atoms in the GB is reminiscent of amorphous packings [34], as evidenced by the large abundance of 555 pairs in the CN analysis, and is also consistent with the results of Swygenhoven et al in that the GB layer is confined to only a few atomic layers.
Two distinct events caused significant decreases in the energy. In the first event, a GB atom obtains local HCP coordination. The HCP core forms a stable structure, which lasts for the remainder of the simulation. The second event is a flattening of the interface regions. The event consists of 20 transitions where 10 atoms leave the GB to obtain FCC coordination and the GB atoms rearrange into larger plateaus. The atomic ordering at the interface between the GB layer and the upper crystal is mostly an FCC(100) plane whereas the interface with the lower crystal is mostly an FCC(111) plane. During the event, most of the activity takes place on the less-stable (001)-like side.
If one had to perform a corresponding simulation with an ordinary MD, this simulation would demand 15 years' exclusive use of a single 2 GHz PowerPC. Our simulation was completed over a period of 10 days where a farm of 100 volunteer processors collaborated, which clearly demonstrates the potential of the AKMC algorithm to handle the timescale problem for certain physical systems.