Prerpints.org logo
Preprint
Article

This version is not peer-reviewed.

Structure and Migration Mechanisms of Small Vacancy Clusters in Cu: A Combined EAM and DFT Study

A peer-reviewed version of this preprint was published in:
Nanomaterials 2023, 13(9), 1464. https://doi.org/10.3390/nano13091464

Submitted:

06 April 2023

Posted:

07 April 2023

You are already at the latest version

Abstract
Voids in face-centered cubic (fcc) metals are commonly assumed to form by the aggregation of vacancies; however, the mechanisms of vacancy clustering and diffusion are not fully understood. In this study, we use computational modeling to provide a detailed insight into the structures and formation energies of primary vacancy clusters, mechanisms and barriers for their migration in bulk copper, and how these properties are affected at simple grain boundaries. The calculations were carried out using Embedded Atom Method (EAM) potentials and Density Functional Theory (DFT) and employed the Site-Occupation Disorder code (SOD), the Activation Relaxation Technique nouveau (ARTn) and the Knowledge Led Master Code (KLMC). We investigate stable structures and migration paths and barriers for clusters of up to six vacancies. Migration of vacancy clusters occurs via hops of individual constituent vacancies with di-vacancies having a significantly smaller migration barrier than mono-vacancies and other clusters. This barrier is further reduced when di-vacancies interact with grain boundaries. This interaction leads to the formation of self-interstitial atoms and introduces significant changes into the boundary structure. Tetra-, penta-, and hexa-vacancy clusters exhibit increasingly complex migration paths and higher barriers than smaller clusters. Finally, the direct comparison with the DFT results shows that EAM can accurately describe the vacancy-induced relaxation effects in the Cu bulk and in grain boundaries. Significant discrepancies between the two methods were found in structures with a higher number of low-coordinated atoms, such as penta-vacancies and di-vacancy absortion by grain boundary. These results will be useful for modeling the mechanisms of diffusion of complex defect structures and provide further insights into the structural evolution of metal films under thermal and mechanical stress.
Keywords: 
vacancies
;  
diffusion
;  
metals
;  
density functional theory
;  
embedded atom method
Subject: 
Chemistry and Materials Science  -   Metals, Alloys and Metallurgy

1. Introduction

Copper is a typical face-centered cubic (fcc) metal with a wide range of applications due to its reliability, low resistivity, and high thermal and electrical conductivity [1,2,3]. It is often used at elevated temperatures or in high-stress conditions [4,5,6]. Voids are commonly formed in fcc metals under high-stress conditions [7,8,9] leading to severe degradation effects [2]. They are thought to be caused by the supersaturation of vacancies (V n , with n being the number of vacancies) [10,11], which diffuse and aggregate either at grain boundaries or triple junctions [3,12,13,14]. Under such conditions, vacancies in fcc metals have been reported to form via two main mechanisms: i) thermally, after quenching, and ii) via the nucleation of two-dimensional defects, such as dislocations. In the first case, quenching takes place at temperatures lower than the melting temperature of the metal [15,16]. In the latter case, migrating dislocations formed via plastic deformation intercept to form steps normal to the slip plane, referred to as jogs. As jogs migrate, rows of vacancies along with lines of interstitial atoms are formed [17,18,19,20,21,22,23]. Experimental studies in Al have shown that clustering of vacancies is already initiated below the room temperature [24].
Since in metals like Cu and Al [12,18,24], such small clusters are not easily identified experimentally, theoretical modeling is often used to determine their characteristics. Theoretical investigations have shown that di-vacancies in fcc metals are significantly more mobile than mono-vacancies or clusters with three or more vacancies [25]. Furthermore, di-vacancy migration has been considered to determine the mobility of larger clusters [26]. Although di-vacancies appear to play an important role in the aggregation of vacancies in fcc metals [27,28], the main mechanism leading to their high mobility in Cu and inside GBs remains unclear. Moreover, theoretical studies of the properties of vacancy clusters in Al indicated that di-vacancies are unstable [29]. When considering the diffusion of larger clusters, tri-vacancies have been reported to be significantly less mobile compared to di-vacancies, serving mostly as nucleation centers [24]. Penta-vacancies have been regarded as stable clusters, with considerably higher binding energies and migration barriers compared to clusters with less than five vacancies [24,30,31]. For hexa- and higher numbers of vacancies, the first perfect stacking fault tetrahedral (SFT) structures have been reported [32].
The role of vacancies in void formation has been highlighted by theoretical simulations. However, the vacancy clustering mechanisms are not fully understood [24]. In particular, grain boundaries (GBs) can serve as efficient sinks for point defects in metals [33,34,35]. Reports have suggested how GBs can reduce the formation, binding, and segregation energies of nearby vacancies [36,37]. Based on experimental and theoretical studies, structural alterations and migration of GBs during recrystallization are considered to be the main contributors to the accumulation of vacancies. Thus, vacancy nucleation and agglomeration can lead to the formation of microstructural defects, such as voids, with other studies describing how strain application can lead to the formation of voids at GBs [38,39,40,41,42]. Despite these insights, the effects of GBs on the diffusion barriers of vacancies remain unclear. In particular, the behavior of di-vacancies inside GBs needs deeper understanding.
One of the main hurdles to elucidating these processes is the need for large-scale simulation models. The investigation of different cluster configurations following different migration paths requires a large number of calculations. Therefore, forcefields (FFs) are often employed [24,25], allowing the simulation of thousands or even millions of atoms. However, depending on the number of parameters and fitting process, FFs simulations can lead to results that deviate from those obtained using ab initio approaches [43,44,45,46,47]. Quantum mechanical simulations, on the other hand, are more accurate but have limitations due to the computational cost. Thus, optimally, a synergy of the two approaches is required [25,48,49].
In this work, we use both the Embedded Atom Method (EAM) and Density Functional Theory (DFT) to identify the most favorable configurations of small vacancy clusters in Cu and study their migration paths and barriers. By combining both methods, we are able to explain why the di-vacancy diffusion barrier is much lower than that for mono-vacancy in fcc Cu. For larger clusters, we use the Site-Occupation Disorder code (SOD) [50] and the Activation Relaxation Technique nouveau (ARTn) [51] to accurately and more reliably identify the mechanisms and barriers for their migration. Since voids and vacancies are strongly related to GBs, the interaction range between vacancies and the [100](210) Σ 5 twin grain boundary is identified and this GB is shown to act as an efficient sink for vacancies and di-vacancies. The most favorable absorption sites along with the absorption region of mono-vacancies in Cu GB are determined. By combining DFT, EAM, and KLMC [52], we demonstrate that the most favorable configurations of di-vacancies inside the GB introduce strong structural changes. Finally, by comparing the results obtained with DFT and EAM we assess the accuracy of EAM in describing the energetic properties of vacancies in Cu. The results show that, although EAM can accurately describe the defect-induced lattice distortions, significant deviations arise in the main energetic properties of the examined clusters. The discrepancy between the two methods stems from the fact that, unlike DFT, EAM simulations favor the formation of self-interstitial atoms (SIAs) in Cu structures containing lower-coordinated sites, such as larger vacancy clusters and GBs with di-vacancies. These results are useful for understanding the behavior of vacancy clusters in fcc metals and provide further insights into the mechanisms of degradation of metal films under thermal and mechanical stress.

2. Methods of Calculations

2.1. Computational Details

2.1.1. Embedded Atom Method

The EAM simulations of vacancy clusters in Cu are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code [53]. For the geometry optimization, the conjugate gradient (CG) optimization algorithm is used with a force tolerance of 10 − 10 eV/ Å. The cells were further equilibrated under an isothermal-isobaric NPT (constant-temperature, constant-pressure) ensemble at constant temperature (0 K) for 0.02 ns using a timestep of 1 fs.
The EAM potential developed by Mishin et al. [54] was used in most calculations. This potential has been fitted semi-empirically and is being widely used for pure Cu systems [55]. Even though many reports highlight the efficiency of such potentials, other studies of fcc metals emphasize the need for developing modified EAM potentials for the accurate description of the properties of point defects and the need for further testing [55,56,57,58].

2.1.2. DFT Calculations

DFT calculations are performed using the CP2K code [59] with the Generalized Gradient Approximation (GGA) Perdew–Burke–Ernzerhof (PBE) functional [60] and at the Γ point. Due to the large size of computational cells, only one k-point is used. CP2K combines the plane wave and Gaussian basis sets (GPW) method [61]. Initially, the cut-off energies for the real-space (RS) integration grid and the mapping of the Gaussians onto the RS grid were converged. An energy cut-off of 550 Ry (7483 eV) is used with a 50 Ry (680 eV) relative cut-off. The Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) was used for geometry optimization [62]. The Fermi–Dirac smearing method is used at a 300 K electronic temperature along with the Broyden mixing method, including 10% of the new density [63] to improve convergence. For Cu atoms, an optimized short-range double zeta for valence electron plus polarization functions (DZVP-MOLOPT-SR-GTH) basis set [64] is used with a Goedecker–Teter–Hutter (GTH) PBE pseudopotential [65].
The simulation cell sizes used must be sufficient to include all the lattice distortion effects induced by vacancy clusters in the process of diffusion. To evaluate the simulation cell size needed in order to include these effects and avoid interactions between periodically translated images of the vacancies, different supercell sizes were tested (2×2×2–6×6×6). It was found that at least a 4×4×4 256-atom supercell was needed for clusters with up to 6 vacancies. Therefore, for the NEB calculations of di-vacancies and tri-vacancies a 4×4×4 256-atom supercell is used. Simulations of clusters with more than three vacancies are carried out in a 5×5×5 500-atom supercell.
To study the interaction of vacancies with GBs, we choose the [100](210) Σ 5 twin boundary, which is one of the lowest energy GBs in Cu [66] and is commonly used in atomistic simulations [67,68,69,70]. The simulation cell was periodically translated in the X, Y and Z. On the Z-axis, an additional 10 Å vacuum was added in order to avoid interactions between periodically translated images. To simulate the absorption of mono- and di- vacancies in the GB, the 304 (15.66 Å×5.56 Å×40.92 Å) and 912 (24.07 Å×14.35 Å×40.92 Å) atom cells are used in DFT and EAM, respectively. To compute the diffusion barriers of di-vacancies towards the GB, two simulation cell sizes are used in DFT, 152 and 304-atoms, whereas five sizes are examined using EAM, namely 152, 304, 456, 608, and 912-atom cells. To render the resulting structures, VESTA [71] and OVITO [72] are used. Grain boundary structures are constructed using Atomsk [73].
Calculations of migration barriers of small Cu vacancy clusters (V 1 , V 2 , V 3 ) through various diffusion paths are performed using the climbing image nudged elastic band (CI-NEB) method implemented in the CP2K code [74]. For all the investigated diffusion paths, five replicas/images are used. Direct inversion in the iterative subspace (DIIS) is used as the optimizer for the band [75].
Formation energies per added vacancy for different clusters are calculated as:
E form = ( E n − E n − 1 ) + μ Cu ,
where E n and E n − 1 are the energies of the crystal with n and n − 1 vacancies. For μ Cu , the energy of a Cu atom in the bulk is used. Binding energies per added vacancy in a cluster are computed as follows:
E bind = ( E n − E n − 1 ) + μ Cu − E 1 f ,
where E 1 f is the formation energy of a mono-vacancy. Positive binding energy values correspond to a favorable addition of a vacancy.
Vacancy absorption energies at GB are computed as the energy difference between a vacancy located in the bulk and at varying positions in the GB. Negative absorption energy values indicate a favorable position for the vacancy.

2.2. Vacancy Cluster Identification

One of the aims of the current study is the identification of the most favorable geometries of vacancy clusters with a different number of vacancies. A cluster is defined as a group of vacancies positioned at a certain distance apart. In order to determine that distance, the interaction range between vacancies needs to be defined. Vacancies have negligible interaction energies (less than 0.05 eV) when separated by more than 5 Å. To reduce the number of configurations of V n clusters by taking into account the symmetry of the system, the Site-Occupation Disorder Code (SOD) is used. The code applies isometric transformations or symmetry operations, like translations, rotations, and reflections, to determine all unique cluster configurations. If an isometric transformation that can convert one cluster configuration into an already identified one is found, the cluster is not considered further [50]. All nonequivalent configurations identified by SOD are first optimized using LAMMPS/EAM, with the lowest energy configurations being reoptimized using CP2K/DFT.
In the bulk, for V 3 , V 4 , V 5 , and V 6 14, 72, 223, and 870 non-equivalent configurations, respectively, are identified using SOD within a 7.23 Å×7.23 Å×7.23 Å cubic 2×2×2 32-atom supercell of fcc Cu. To determine the structure and energy of each cluster, 500-atom cells containing the cluster are fully relaxed using the EAM potential, and the lowest energy configurations are identified within the range of 10–25 kT, depending on the number of vacancies, where k is the Boltzmann’s constant and T is the room temperature. We consider this energy range sufficient to take into account the thermodynamically accessible structures. The geometries of these structures were further optimized using DFT to identify the lowest energy configurations.
In regions with low symmetry, such as GBs, the Knowledge Led Master Code (KLMC) is employed. KLMC is a structure prediction method enabling us to refine input files and generate outputs that can be utilized by third-party electronic structure software [52,76,77]. A detailed description of how the code operates can be found in ref. [52]. Using the code, one can predict structures not only in the bulk but also at free surfaces. In this study, KLMC is used to identify all possible di-vacancy configurations in the GB. Initially, 1500 V 2 configurations generated via KLMC are optimized with EAM. The lowest energy structures within the range of 10 kT are reoptimized using DFT.

2.3. Calculation of Diffusion Barriers of Vacancy Clusters

The energetically most favorable cluster configurations in the bulk are used as the initial and final configurations for the DFT diffusion barrier calculations of mono-, di- and tri-vacancies using the nudged elastic band CI-NEB method [74]. For such clusters, due to the crystal symmetry, CI-NEB is sufficient to identify all the different diffusion paths. The path with the lowest diffusion barrier is defined as the most favorable.
For larger clusters, the saddle points cannot be identified based on crystal symmetry. To identify and compute the migration barriers of V 4 , and V 5 clusters, the Activation Relaxation Technique nouveau [51,78] is used in tandem with EAM and DFT. In ARTn, nearby local minima are searched by randomly displacing atoms within a certain radius (deformation radius) by a predetermined distance. To identify the cluster migration paths along with the transition states, a three-step process is followed. At the start of the activation, atoms within the deformation radius with respect to the atom closest to the center of mass of the vacancy cluster of interest are displaced along random directions. To achieve that, a nonzero term in the force is created in order for the configuration to escape the harmonic well. The configuration is then shifted by iteratively applying a redefined force to a neighboring saddle point and converged to a new local minimum [79]. For the new local minimum point to be accepted, the energy difference between the two local minima needs to be lower than a certain energy threshold. The same process is repeated for the identification of more saddle points [78,79,80] until the system reaches an equivalent cluster configuration with the center of mass of the cluster displaced compared to the initial configuration. This three-step process is repeated several times to identify different migration paths. A migration path typically includes several individual vacancy moves and corresponding barriers. Once all the paths are identified, we define the most favorable path as the one with the lowest barriers.
More specifically, each search starts from the lowest energy configuration of each cluster. When the first saddle point is found, the structure is fully optimized using EAM. After the optimization, a second local minimum point is identified using the ARTn code. The search is repeated until the nearest equivalent cluster configuration is found. All the minima and saddle point configurations identified via ARTn/EAM for the lowest energy path are reoptimized using DFT. For the saddle point DFT optimization, the first nearest neighbor atoms from the introduced vacancies are fixed. In the case of V 6 , since more than 20 configurations are identified using ARTn, only EAM was used for their optimization. To define the accuracy of the saddle point identification process, the deformation radius is chosen at 5 Å. A criterion to accept or reject new minimum points, also referred to as events, is chosen at 0.5 eV. The maximum distance for breaking the crystal symmetry by random displacement of atoms within the deformation radius is set at 0.05 Å.
For V 4 , five runs are conducted, whereas for larger clusters (V 5 , V 6 ) the number of runs is increased to ten. For each run, a maximum number of one hundred events are generated, with acceptance event rates ranging between 91–93%. Following this process allows us to understand the initial steps for the deformation and migration of vacancy clusters and preferable migration paths.

3. Results of Calculations

The same color code is used in all the figures below, with blue atoms corresponding to Cu atoms, vacancies shown in red, and SIAs are shown in white.

3.1. Structure of Vacancy Cluster

Figure 1 (a) and (b) illustrate the relative and formation/binding energies, respectively, of the lowest energy configuration of V 2 - V 6 clusters identified using the SOD code and optimized with DFT. For the formation and binding energies in Figure 1 (b), the EAM results are also included. Figure 1 (c) shows the lowest energy configurations of V 3 - V 6 clusters optimized using DFT. The configurations are shown in the same energetic order as in Figure 1 (a). As expected, vacancies tend to cluster in closed-packed structures. EAM predicted the same lowest-energy tri-vacancy configuration as DFT, which forms the closest packed cluster with all three vacancies relaxing at first nearest neighbor (1NN) distances from each other, as shown in Figure 1 (c) V 3 5. The difference in energies between the lowest and less compact second-lowest tri-vacancy clusters is approximately 0.1 eV.
EAM and DFT identified the same lowest energy configuration for a tetra-vacancy, which is the close-packed cluster shown in Figure 1 (c) V 4 5. The energy difference between the two lowest energy V 4 structures seen in Figure 1 (a) optimized with DFT is less than 0.1 eV. These results indicate that rapid transitions between different configurations of V 3 and V 4 clusters should be taking place even at relatively low temperatures.
The lowest energy configuration of the V 5 cluster identified using DFT is shown in Figure 1 (c), where V 5 5 forms the closest packed penta-vacancy. We note that this is the only case where EAM predicted a different lowest energy configuration compared to DFT. According to EAM, the lowest energy penta-vacancy cluster consists of six vacancies and a SIA located in the geometric center of the cluster, as shown in Figure 1 (c) V 5 1. A similar structure has been previously reported in theoretical studies of Al and Cu using other forcefields [24]. Our EAM calculations predict that this structure is lower in energy compared to the second-lowest energy structure by 0.13 eV. However, according to DFT, this configuration is 0.67 eV higher in energy than the pyramidal V 5 5 configuration.
Both EAM and DFT predicted the same lowest energy V 6 cluster configuration. As seen in Figure 1 (c) V 6 5, this configuration forms a symmetric octahedral void, the same as the V 5 configuration shown in Figure 1 (c) V 5 1 but without the SIA. Both the V 5 and V 6 configurations shown in Figure 1 (c) V 5 1 and V 6 4, respectively, are the only ones including SIAs predicted by DFT. In the case of V 6 , this cluster using DFT is 0.14 eV higher in energy compared to the lowest V 6 configuration. It consists of three SIAs and nine vacancies. Similar lowest energy V 6 configurations have been reported in previous DFT studies [81]. These results indicate that the symmetry of the cluster plays an important role in the stabilization of SIAs.
We note that previous studies that employed forcefields have emphasized the role of SIAs in the microstructural evolution of bcc metals, like Fe and W, due to their high mobility [82,83]. Our results demonstrate that EAM more readily predicts the formation of SIAs than DFT. A number of stable configurations of V 4 , V 5 , and V 6 clusters including SIAs obtained with EAM were found to be unstable with DFT. This suggests that EAM potentials are not always accurate for predicting the behavior of SIAs.

3.2. Formation and Binding Energies of Vacancies

Previous MD studies [18] illustrated that, as jogs nucleate and dislocations glide, the formed rows of vacancies also include discontinuities. The formation of a vacancy cluster along the path that jogs migrate may affect the distributions of individual vacancies. Thus, it is important to determine the range and strength of the interaction between mono-vacancies and vacancy clusters.
The DFT and EAM results, summarized in Figure 1 (b), show similar trends for binding and formation energies per added vacancy. These are defined by equations 1 and 2 above. The plot includes the binding/formation energies of the lowest energy configurations for each cluster obtained through DFT, as seen in Figure 1 (c). As the number of added vacancies increases, the binding and formation energies per vacancy increase and decrease, respectively. The penta-vacancy formation energy is slightly higher per vacancy than that for the tetra-vacancy due to the less compact configuration. The addition of the sixth vacancy creates a compact and symmetric configuration, significantly reducing the formation energy per vacancy. We need to point out that our DFT and EAM results predict the formation of a di-vacancy as energetically favorable, unlike previous studies in fcc metals [29,84]. These results indicate that, as the number of attached vacancies in small clusters increases, the addition of the following vacancy becomes more favorable than the decomposition of the cluster. It is still not clear up to what number of vacancies the addition of further vacancies would start being unfavorable. However, such an investigation is beyond the scope of the current study.

3.3. Cluster Migration Paths in the Bulk Cu

3.3.1. Migration of V 1 – V 3 clusters

Both EAM and DFT calculations demonstrate that individual mono-vacancies have a very weak (lower than 0.1 eV) interaction with other vacancies and small clusters V 1 –V 5 at distances between the centers of symmetry of the clusters and the mono-vacancy exceeding 5 Å. The most significant interaction takes place at distances below 3.6 Å, indicating that these interactions are very short-ranged. Hence, diffusion of clusters and vacancies is the only way for them to meet in perfect lattice to form bigger clusters and voids. To investigate the migration mechanisms of small clusters, different migration paths have been calculated for each cluster using NEB, but only the lowest energy barriers are included in Figure 2 (a).
The mono-vacancy prefers to hop between the 1NN sites with a barrier of 0.65 eV, in good agreement with previous computational studies reporting barriers of 0.69 eV [85]. More paths are available for a di-vacancy. The paths that required the dissociation of the V 2 cluster into two individual mono-vacancies at a next nearest neighbor (NNN) distance are found to have a barrier of 0.55 eV or above. However, when one of the mono-vacancies hops to an adjacent site remaining at the 1NN distance from the other vacancy, the barrier is much smaller at 0.4 eV. The main reason for this reduction is not only due to the non-dissociation of the cluster but also to the reduced crystal distortion energy. When an individual vacancy migrates, an adjacent Cu atom occupies the vacancy site (Figure 2 (b)). This atom has to diffuse through a gate formed by four Cu atoms for such a move to take place. During this process, each of these four Cu atoms is displaced by 0.2 Å. When the two vacancies remain at the 1NN distance during the migration, the number of atoms that need to be displaced is reduced by one. Thus, the energy barrier for cluster diffusion is reduced. This explains why in fcc metals di-vacancies are considered to be more mobile than other clusters.
The migration mechanism for V 3 shown in Figure 2 (c) involves two steps: first, one of the vacancies hops into a nearby site, forming a less compact transient configuration of three vacancies. Secondly, another vacancy hops to restore the compact rectangular configuration at a new position of the center of mass (see Figure 2 (c)). This mechanism involves only single vacancy hops and has a barrier of 0.52 eV.

3.3.2. Elementary migration paths of V 4 –V 6 clusters

In small clusters, the single vacancy hop mechanisms resulted in lower barriers compared to a collective motion of the cluster. In this section, we consider whether this mechanism also works for larger clusters. Since the diffusion of larger clusters involves a substantially larger number of steps, the ARTn code is used to search for adjacent local minima and the corresponding energy barriers. The search continues until a migration path for transitioning from the minimum energy structure found using SOD to the same structure at a different nearby site in the crystal is identified. As the acceptance criterion for new local minima in ARTn we used 0.5 eV and only single vacancy hopping mechanisms were identified.
Barriers identified for migration of the V 4 cluster and some of the intermediate configurations are shown in Figure 3 (a) and (b) (i), respectively. This mechanism involves consecutive jumps of single vacancies, with a DFT calculated barrier of 0.84 eV. As can be seen in Figure 3, this involves a three-step process between two almost equivalent minima which proceeds through two transient configurations and three saddle points. The transient configurations during the migration can be seen in Figure 3 (b). V 4 migration involves the same diffusion mechanism as the one reported for V 2 (see Figure 2 (b)). The computed barriers are close to the results of previous studies for fcc Ni, where the 0.66 eV barrier was calculated for tetra-vacancies using EAM [25].
Figure 3 also shows the migration path and barriers for the lowest energy V 5 cluster. For identifying the diffusion path and prior to optimizing with DFT, the ARTn search employed the EAM potential developed by Ackland et al. [86]. This potential is used since it predicted the same lowest energy penta-vacancy configuration as DFT (see Figure 1 (c) V 5 5). The highest barrier for this cluster is approximately 0.84 eV. Similarly high barriers were found for V 5 in Al [87] using EAM, namely 1.22 eV. Therefore, such clusters will be mobile only at high temperatures. The migration process involved nine configurations: two global minimum points, four local minimum points, and four saddle points. Each step involved the migration of one vacancy (see Figure 3 (b) (ii)).
The migration path for the V 6 cluster is more complicated. Using ARTn and EAM we have identified 18 transition configurations with the highest barrier of 0.96 eV. The calculated barriers for all clusters are summarized in Table 1. The highest barrier for the most preferred path is regarded as the diffusion barrier of the cluster. They indicate that migration of V 4 , V 5 , and V 6 clusters may become important only at high temperatures reaching 1000 K, as has recently been demonstrated for Ni [88].

3.3.3. Crowdion motion of vacancy clusters

Previous studies have mentioned the collective one-dimensional (1D) crowdion motion as a favorable diffusion mechanism for vacancy clusters in metals, especially in bcc lattices [89,90,91]. Our DFT results show that the barrier for crowdion motion of a di-vacancy toward the Cu(111) surface is higher by 0.16 eV than the identified lowest barrier of a V 2 . The most favorable 1D configuration of a V 3 cluster is higher in energy compared to the lowest energy V 3 configuration by approximately 0.2 eV and 0.15 eV with DFT and EAM, respectively. The crowdion displacement barrier for V 3 cluster is by 0.15 eV higher than the most favorable V 3 barrier.
EAM also showed that linear crowdion configurations of V 4 , V 5 , and V 6 clusters were higher in energy by 0.46 eV, 0.78 eV, and 1.22 eV, respectively, compared to the lowest energy V 4 , V 5 , and V 6 clusters. The same energetic differences using DFT were 0.2, 0.55, and 0.92 eV for V 4 , V 5 , and V 6 , respectively. These results indicate that within both EAM and DFT, V 4 , V 5 , and V 6 prefer to form close-packed clusters. Therefore, for clusters with a small number of vacancies, a 1D crowdion diffusion would be unfavorable. However, since a series of vacancies and SIAs are expected to form as jog-connected dislocations migrate, our results do not rule out the possibility of crowdion motion of a 1D linear configuration of clusters of more than six vacancies.

4. Vacancy Interaction with Grain Boundaries

Our calculations of vacancy migration paths in the bulk Cu have demonstrated that the migration of larger clusters (V 4 - V 6 ) will require substantially higher energy barriers compared to smaller clusters. However, under thermodynamic equilibrium conditions, such clusters are unlikely to form in perfect Cu due to the high vacancy formation energies [92] even at high temperatures [48]. We investigate how the properties of vacancies are affected under non-equilibrium conditions where 2D defects are present. In particular, vacancies and vacancy clusters can diffuse to GBs and become absorbed by them with a varying degree of efficiency.

4.1. Interaction of mono-vacancy with grain boundary

Figure 4 (a) (i) shows the Σ 5 twin grain boundary used to determine how the presence of two-dimensional defects that act as efficient sinks can affect the mobility of vacancies. By introducing a vacancy in each of the numbered sites, the most favorable absorption sites for a vacancy in the GB are identified. The selection of sites was based on previous DFT segregation studies in Cu GB [93,94]. The simulation cells for DFT and EAM comprised of 304 and 912 atoms, respectively, are shown in Figure 4 (a) (ii). EAM and DFT results are included in the graph shown in Figure 4 (b). The plot includes the computed absorption energies both prior to (single point) and after the geometry relaxation at various sites in the GB. Prior to relaxation, EAM and DFT are in good agreement in identifying the equivalent sites 5 and 7 as the most favorable positions, whereas site 6 is found to be unfavorable. The same most favorable segregation sites have been identified in ref. [68] using DFT. We note that our 304-atom cells are considerably larger compared to 76–112-atom cells used in previous similar DFT studies [95,96,97], which allows us to carry out full relaxation of the GB structure as a result of the absorption of vacancies. After the relaxation, both methods show that the vacancy migration from sites 3/9, 4/8 to site 5/7 is accompanied by an energy gain of around 1.0 eV. DFT predicted a stronger GB–vacancy interaction by approximately 0.1 eV compared to EAM.
The interaction range between vacancies and GBs (absorption region) is demonstrated in Figure 4 (c), where the absorption energies of a mono-vacancy at different distances from the center of symmetry of the GB are calculated. The absorption region is approximately 6.9 Å wide, corresponding to an interaction radius of 3.45 Å between vacancies and GBs. The similar absorption region of 6.36 Å was identified in ref. [39] for mono-vacancies in Σ 5 grain boundary in bcc W using EAM. DFT and EAM absorption energies are in good agreement, with both methods predicting a gain of 1.0 eV for the absorption of a mono-vacancy. These results are broadly in agreement with the calculations presented in ref. [98].

4.2. Interaction of di-vacancy with grain boundary

Our results, as well as previous data [25], predict the higher mobility of di-vacancies than mono-vacancies in the bulk of Cu. Hence, it is of interest to study the interaction of di-vacancies with the regions that are prone to absorb point defects, such as GBs. First, we search for favorable absorption sites for di-vacancies within 6.9 Å from the GB. Different configurations were tried using a similar process to the one used for the vacancy cluster investigation in the bulk. Due to the larger number of lower-coordinated atoms in GB compared to the bulk, KLMC was used to identify all different configurations of di-vacancies. All identified V 2 configurations were optimized with EAM, with the lowest energy structures being reoptimized with DFT.
Figure 5 (a) shows the lowest energy initial V 2 configurations prior to relaxation using DFT. The configurations are shown from (i) to (v) in energetic order, with (v) resulting in the lowest energy GB configuration after the relaxation. Both DFT and EAM predict the initial configurations of the two vacancies sitting at a 1NN distance 2.6 Å apart, to be unfavorable. DFT identifies two vacancies initially separated by 3.7 Å as the most favorable initial V 2 configuration (see Figure 5 (b) (v)). Hence, di-vacancies tend to dissociate close to this GB. Such behavior has previously only been identified for larger vacancy clusters where nanovoids were found to dissociate when interacting with migrating GBs [34].
Figure 5 (b) and (c) include the lowest energy fully relaxed GB configuration with a di-vacancy using EAM and DFT, respectively. The incorporation of a di-vacancy into the GB causes significant relaxation effects that facilitate the formation of SIAs. In the case of EAM, (Figure 5 (b)), the geometry optimization of the periodically translated GB with di-vacancy results in the shift of the main symmetry axis of the GB. Also, Σ 5 deltoids are transformed into filled deltoids. Similar structural changes have been reported in ref. [99]. On the other hand, in the DFT relaxed GB configuration with a di-vacancy, seen in Figure 5 (c), the initial Σ 5 deltoids are transformed into a split-deltoid configuration through the formation of SIAs that sit within the excess free volume in the GB area. We note that significant differences in the relaxed GB structures are seen in Figure 5 (b) and (c). Like in the bulk, a considerably higher number of SIAs are formed within the GBs when using EAM (Figure 5 (b)) compared to DFT (Figure 5 (c)). The latter illustrates that EAM could not capture the same relaxation effects induced by the di-vacancy when compared with DFT. This stems from the noted above trend that EAM simulations favor the SIAs formation in Cu structures containing lower-coordinated sites, such as larger vacancy clusters and GBs with di-vacancies.
Finally, we employed NEB calculations to investigate the di-vacancy absorption into the GB. Figure 6 shows the V 2 migration path and energy profile obtained using NEB with five images. The initial configuration corresponds to a di-vacancy in the bulk but within the absorption region of 6.9 Å defined in Figure 4 (c). The final configuration is the lowest energy configuration of a di-vacancy inside the GB seen in Figure 5 (a) (v). The di-vacancy configurations during the absorption process are shown in Figure 6 (b) in the same order as in Figure 6 (a). DFT and EAM predict the lower than 0.05 eV difference between configurations corresponding to images 1 and 2. The absorption of the di-vacancy by GB results in the 2.0 eV gain in energy caused by the strong structural relaxation. The high 2.0 eV reverse barrier illustrates the strong binding energy of di-vacancy to GB.
This strong lattice relaxation prompted us to test how the size of the periodic cell affects the amount of energy gained when a di-vacancy moves toward the GB. Two simulation cell sizes are tested with DFT and four with EAM. The largest GB cell used for DFT contains 304 atoms, whereas, for EAM, the largest cell contains 606 atoms. As seen in Figure 6 (a), using both DFT and EAM, 304-atom or larger cells are required to achieve converged energy values. However, for this cell size, the barriers computed with both methods are in good agreement.

5. Conclusions

We have investigated the structure, migration paths and barriers of vacancy clusters in the bulk of fcc Cu and near the Σ 5 twin grain boundary. In the bulk, several most stable cluster structures up to six vacancies have been identified. The most significant interaction between vacancies takes place within one lattice parameter distance. As expected, the DFT results in the bulk show that formation and cluster binding energies per vacancy reduce and increase, respectively, and the addition of a further vacancy becomes more favorable than the decomposition of the cluster as the number of vacancies increases.
Using both DFT and EAM, the migration paths and barriers for clusters of up to three vacancies were computed. To identify the elementary diffusion mechanism and calculate the migration barriers of larger clusters, DFT and EAM were used with ARTn. In the bulk, di-vacancies have the lowest migration barrier amongst the examined clusters, which confirms previous studies in fcc metals. V 4 migration was found to follow the same diffusion mechanism as that of di-vacancies. V 4 , V 5 , and V 6 clusters have considerably higher migration barriers compared to V 2 .
The diffusion barrier of a di-vacancy was further reduced in the vicinity of the Σ 5 twin grain boundary. The absorption of V 2 by the GB induced the formation of SIAs within the free volume formed by the GB. In addition, the GB structure was found to undergo significant changes from the initial deltoid-like structure due to the incorporation of di-vacancy. Our results also showed that di-vacancies prefer to dissociate inside the GB, a behavior that was previously shown only for larger vacancy clusters. Hence, these results indicate that for the initiation of voids in these regions, the presence of other factors may need to be considered, such as additional point defects or impurities.
Finally, since a number of studies in the past employed forcefield methods to model of vacancy clusters in fcc metals, we evaluated the accuracy of EAM compared to DFT calculations. Both EAM and DFT demonstrate that clusters prefer to form closed-packed structures. Excluding V 5 , the same lowest energy cluster configurations were identified using DFT and EAM. The EAM potential [54] proved to be reliable for describing the binding and formation energies of vacancies in the bulk of Cu. The EAM results are in agreement with DFT regarding the lowest energy structures of V 3 , V 4 , and V 6 clusters within the bulk and the description of the behavior of vacancies within the GB and the vacancy induced relaxation effects. Also, EAM results are in good agreement with DFT in predicting the energy gain due to the migration of a di-vacancy in the GB, provided a sufficiently large periodic cell is employed. However, when compared with DFT, EAM did not capture the relaxation effects induced in the GB due to the presence of a di-vacancy. More significant discrepancies between the DFT and EAM results were observed in cases where structures with a high number of low-coordinated sites were considered, such as in penta-vacancies and grain boundaries.
These results are not only useful for our understanding the behavior of vacancy clusters in fcc metals, but also provide a direct comparison between EAM and DFT methods that can be used as a reference for future studies. Our approach can be employed to provide a better understanding of structure and diffusion properties of more complex structures, such as complexes of vacancies with metallic and non-metallic impurities.

Author Contributions

Conceptualization, V.F., M.K., E.K. and A.L.S; methodology, V.F. and D.M.-F.; software, V.F. and D.M.-F.; validation, V.F., D.M.-F., R.B. and L.R.; formal analysis, V.F., D.M.-F. and A.L.S.; investigation, V.F., D.M.-F and A.L.S.; resources, V.F. and A.L.S.; data curation, V.F.; writing—original draft preparation, V.F., D.M.-F. and A.L.S.; writing—review and editing, V.F., D.M.-F., M.K., R.B., L.R., E.K. and A.L.S.; visualization, V.F. and R.B.; supervision, A.L.S.; project administration, V.F. and A.L.S.; funding acquisition, V.F. and A.L.S. All authors have read and agreed to the published version of the manuscript.

Funding

D.M.-F. and A.L.S. acknowledge funding by EPSRC (grant EP/P013503/1). V.F. would like to acknowledge funding by EPSRC (grant EP/L015862/1) as part of the CDT in molecular modeling and materials science. Computational resources on Archer (http://www.archer.ac.uk) were provided via our membership of the UK’s HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202, EP/R029431). The work by E.K. and M.K. was partially funded by the Austrian Research Promotion Agency (FFG, Project No. 881110). R.B. and L.R. gratefully acknowledge the financial support under the scope of the COMET program within the K2 Center “Integrated Computational Material, Process and Product Engineering (IC-MPPE)” (Project No. 859480). This program was supported by the Austrian Federal Ministries for Climate Action, Environment, Energy, Mobility, Innovation and Technology (BMK) and for Digital and Economic Affairs (BMDW), represented by the Austrian research funding association (FFG), and the federal states of Styria, Upper Austria, and Tyrol. L.R. also acknowledges financial support by the Austrian Science Fund (FWF) (P 34179-N). The work by E.K. and M.K. was partially funded by the Austrian Research Promotion Agency (FFG, Project No. 881110).

Institutional Review Board Statement

Not applicable.

Acknowledgments

V.F. and A.L.S. would like to thank Jack Strand, Tom Durrant and David Schmidt for useful comments and help in calculations.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
DFT Density Functional Theory
EAM Embedded Atom Method
LAMMPS Atomic/Molecular Massively Parallel Simulator
ARTn Activation Relaxation Technique nouveau
KLMC Knowledge Led Master Code
PBE Perdew–Burke–Ernzerhof
GPW Gaussian Planewave
GTH Goedecker–Teter–Hutter
NEB nudged elastic band
CI-NEB climbing image nudged elastic band
GB Grain Boundary
SFT Stacking Fault Tetrahedral
fcc face-centered cubic
FF forcefield
NNN Next Nearest Neighbor
1NN First Nearest Neighbor
SIA Self Interstitial Atom
ASE Atomic Simulation Environment
CG conjugate gradient
GGA Generalized Gradient Approximation
RS real-space
BFGS Broyden–Fletcher–Goldfarb–Shanno
DIIS direct inversion in the iterative subspace
E f o r m formation energy
E b i n d binding energy
E 1 f mono-vacancy formation energy
E 1 f formation energy of a cluster with n vacancies
V n vacancy cluster with n vacancies

References

  1. Needleman, A. A Continuum Model for Void Nucleation by Inclusion Debonding. Journal of Applied Mechanics 1987, 54. [Google Scholar] [CrossRef]
  2. Feng, H.; Cheng, M.; Wang, Y.; Chang, S.; Wang, Y.; Wan, C. Mechanism for Cu void defect on various electroplated film conditions. Thin Solid Films 2006, 498, 56–59. [Google Scholar] [CrossRef]
  3. Gondcharton, P.; Imbert, B.; Benaissa, L.; Verdier, M. Voiding phenomena in copper-copper bonded structures: role of creep. ECS Journal of Solid State Science and Technology 2015, 4, P77. [Google Scholar] [CrossRef]
  4. Sharif, A. All-Copper Interconnects for High-Temperature Applications. Harsh Environment Electronics: Interconnect Materials and Performance Assessment 2019.
  5. Lee, T.K.; Chen, Z.; Baty, G.; Bieler, T.R.; Kim, C.U. Impact of an Elevated Temperature Environment on Sn-Ag-Cu Interconnect Board Level High-G Mechanical Shock Performance. Journal of Electronic Materials 2016, 45, 6177–6183. [Google Scholar] [CrossRef]
  6. Li, Z.; Tian, Y.; Teng, C.; Cao, H. Recent Advances in Barrier Layer of Cu Interconnects. Materials 2020, 13, 5049. [Google Scholar] [CrossRef] [PubMed]
  7. Gerstein, G.; Besserer, H.B.; Nürnberger, F.; Barrales-Mora, L.A.; Shvindlerman, L.S.; Estrin, Y.; Maier, H.J. Formation and growth of voids in dual-phase steel at microscale and nanoscale levels. Journal of Materials Science 2017, 52, 4234–4243. [Google Scholar] [CrossRef]
  8. Pang, W.W.; Zhang, P.; Zhang, G.C.; Xu, A.G.; Zhao, X.G. Dislocation creation and void nucleation in FCC ductile metals under tensile loading: A general microscopic picture. Scientific reports 2014, 4, 1–7. [Google Scholar] [CrossRef]
  9. Qi, Y.; Chen, X.; Feng, M. Effect of void defect on c-axis deformation of single-crystal Ti under uniaxial stress conditions: Evolution of tension twinning and dislocations. Journal of Materials Research 2019, 34, 3699–3706. [Google Scholar] [CrossRef]
  10. Katz, J.L.; Wiedersich, H. Nucleation of voids in materials supersaturated with vacancies and interstitials. The Journal of Chemical Physics 1971, 55, 1414–1425. [Google Scholar] [CrossRef]
  11. Fischer, F.; Svoboda, J. Void growth due to vacancy supersaturation–A non-equilibrium thermodynamics study. Scripta Materialia 2008, 58, 93–95. [Google Scholar] [CrossRef]
  12. Kovács, Z.; Chinh, N.Q. Up-hill diffusion of solute atoms towards slipped grain boundaries: A possible reason of decomposition due to severe plastic deformation. Scripta Materialia 2020, 188, 285–289. [Google Scholar] [CrossRef]
  13. Adlakha, I.; Solanki, K. Structural stability and energetics of grain boundary triple junctions in face centered cubic materials. Scientific reports 2015, 5, 1–7. [Google Scholar] [CrossRef] [PubMed]
  14. Adlakha, I.; Solanki, K. Atomic-scale investigation of triple junction role on defects binding energetics and structural stability in α-Fe. Acta Materialia 2016, 118, 64–76. [Google Scholar] [CrossRef]
  15. Liu, Y.N.; Ahlgren, T.; Bukonte, L.; Nordlund, K.; Shu, X.; Yu, Y.; Li, X.C.; Lu, G.H. Mechanism of vacancy formation induced by hydrogen in tungsten. AIP advances 2013, 3, 122111. [Google Scholar] [CrossRef]
  16. Rasch, K.; Siegel, R.; Schultz, H. Quenching and recovery investigations of vacancies in tungsten. Philosophical Magazine A 1980, 41, 91–117. [Google Scholar] [CrossRef]
  17. Yuan, S.; Huang, M.; Zhu, Y.; Li, Z. A dislocation climb/glide coupled crystal plasticity constitutive model and its finite element implementation. Mechanics of Materials 2018, 118, 44–61. [Google Scholar] [CrossRef]
  18. Xu, D.;Wang, H.; Yang, R.; Veyssière, P. Point defect formation by dislocation reactions in TiAl. In Proceedings of the IOP Conference Series: Materials Science and Engineering. IOP Publishing, 2009, Vol. 3, p. 012024.
  19. Smoluchowski, R. Dislocations in ionic crystals (structure, charge effects and interaction with impurities). Le Journal de Physique Colloques 1966, 27, C3–3. [Google Scholar] [CrossRef]
  20. Niu, X.; Luo, T.; Lu, J.; Xiang, Y. Dislocation climb models from atomistic scheme to dislocation dynamics. Journal of the Mechanics and Physics of Solids 2017, 99, 242–258. [Google Scholar] [CrossRef]
  21. Zhou, S.; Preston, D.; Lomdahl, P.; Beazley, D. Large-scale molecular dynamics simulations of dislocation intersection in copper. Science 1998, 279, 1525–1527. [Google Scholar] [CrossRef]
  22. Srivastava, K.; Rao, S.I.; El-Awady, J.A. Unveiling the role of super-jogs and dislocation induced atomic-shuffling on controlling plasticity in magnesium. Acta Materialia 2018, 161, 182–193. [Google Scholar] [CrossRef]
  23. Dupraz, M.; Sun, Z.; Brandl, C.; Van Swygenhoven, H. Dislocation interactions at reduced strain rates in atomistic simulations of nanocrystalline Al. Acta Materialia 2018, 144, 68–79. [Google Scholar] [CrossRef]
  24. Wang, H.; Rodney, D.; Xu, D.; Yang, R.; Veyssière, P. Defect kinetics on experimental timescales using atomistic simulations. Philosophical Magazine 2013, 93, 186–202. [Google Scholar] [CrossRef]
  25. Mahmoud, S.; Trochet, M.; Restrepo, O.A.; Mousseau, N. Study of point defects diffusion in nickel using kinetic activation-relaxation technique. Acta Materialia 2018, 144, 679–690. [Google Scholar] [CrossRef]
  26. Kurishita, H.; Amano, Y.; Kobayashi, S.; Nakai, K.; Arakawa, H.; Hiraoka, Y.; Takida, T.; Takebe, K.; Matsui, H. Development of ultra-fine grained W–TiC and their mechanical properties for fusion applications. Journal of nuclear Materials 2007, 367, 1453–1457. [Google Scholar] [CrossRef]
  27. McLellan, R.; Angel, Y. The thermodynamics of vacancy formation in fcc metals. Acta metallurgica et materialia 1995, 43, 3721–3725. [Google Scholar] [CrossRef]
  28. Bartdorff, D.; Neumann, G.; Reimers, P. Self-diffusion of 64Cu in copper single crystals Monovacancy and divacancy contributions. Philosophical magazine A 1978, 38, 157–165. [Google Scholar] [CrossRef]
  29. Ho, G.; Ong, M.T.; Caspersen, K.J.; Carter, E.A. Energetics and kinetics of vacancy diffusion and aggregation in shocked aluminium via orbital-free density functional theory. Physical Chemistry Chemical Physics 2007, 9, 4951–4966. [Google Scholar] [CrossRef] [PubMed]
  30. Vineyard, G.H. General introduction. Discussions of the Faraday Society 1961, 31, 7–23. [Google Scholar] [CrossRef]
  31. Shimomura, Y.; Nishiguchi, R. Vacancy clustering to faulted loop, stacking fault tetrahedron and void in fcc metals. Radiation Effects and Defects in Solids 1997, 141, 311–324. [Google Scholar] [CrossRef]
  32. Xv, H.; Zhao, J.; Ye, F.; Tong, K. Strain-induced transformation between vacancy voids and stacking fault tetrahedra in Cu. Computational Materials Science 2019, 158, 359–368. [Google Scholar] [CrossRef]
  33. Zhang, L.; Lu, C.; Pei, L.; Zhao, X.; Zhang, J.; Tieu, K. Evaluation of mechanical properties of Σ5 (210)/[001] tilt grain boundary with self-interstitial atoms by molecular dynamics simulation. Journal of Nanomaterials 2017, 12, 1–11. [Google Scholar] [CrossRef]
  34. Zhang, L.; Shibuta, Y.; Lu, C.; Huang, X. Interaction between nano-voids and migrating grain boundary by molecular dynamics simulation. Acta Materialia 2019, 173, 206–224. [Google Scholar] [CrossRef]
  35. Chimi, Y.; Iwase, A.; Ishikawa, N.; Kobiyama, M.; Inami, T.; Okuda, S. Accumulation and recovery of defects in ion-irradiated nanocrystalline gold. Journal of Nuclear Materials 2001, 297, 355–357. [Google Scholar] [CrossRef]
  36. Mehrer, H.; Seeger, A. Interpretation of Self-Diffusion and Vacancy Properties in Copper. physica status solidi (b) 1969, 35, 313–328. [Google Scholar] [CrossRef]
  37. Balogh, Z.; Schmitz, G. Diffusion in metals and alloys. In Physical Metallurgy; Elsevier, 2014; pp. 387–559.
  38. Wan, L.; Geng, W.T.; Ishii, A.; Du, J.P.; Mei, Q.; Ishikawa, N.; Kimizuka, H.; Ogata, S. Hydrogen embrittlement controlled by reaction of dislocation with grain boundary in alpha-iron. International Journal of Plasticity 2019, 112, 206–219. [Google Scholar] [CrossRef]
  39. Chen, N.; Niu, L.L.; Zhang, Y.; Shu, X.; Zhou, H.B.; Jin, S.; Ran, G.; Lu, G.H.; Gao, F. Energetics of vacancy segregation to [100] symmetric tilt grain boundaries in bcc tungsten. Scientific reports 2016, 6, 1–12. [Google Scholar] [CrossRef]
  40. Xu, W.; Zhang, Y.; Cheng, G.; Jian, W.; Millett, P.C.; Koch, C.C.; Mathaudhu, S.N.; Zhu, Y. In-situ atomic-scale observation of irradiation-induced void formation. Nature communications 2013, 4, 1–6. [Google Scholar] [CrossRef]
  41. Uberuaga, B.; Hoagland, R.; Voter, A.; Valone, S. Direct transformation of vacancy voids to stacking fault tetrahedra. Physical review letters 2007, 99, 135501. [Google Scholar] [CrossRef]
  42. Granberg, F.; Nordlund, K.; Ullah, M.W.; Jin, K.; Lu, C.; Bei, H.; Wang, L.; Djurabekova, F.; Weber, W.; Zhang, Y. Mechanism of radiation damage reduction in equiatomic multicomponent single phase alloys. Physical review letters 2016, 116, 135504. [Google Scholar] [CrossRef]
  43. Fröhlking, T.; Bernetti, M.; Calonaci, N.; Bussi, G. Toward empirical force fields that match experimental observables. The Journal of chemical physics 2020, 152, 230902. [Google Scholar] [CrossRef]
  44. Kashefolgheta, S.; Oliveira, M.P.; Rieder, S.R.; Horta, B.A.; Acree Jr, W.E.; Hünenberger, P.H. Evaluating Classical Force Fields against Experimental Cross-Solvation Free Energies. Journal of Chemical Theory and Computation 2020, 16, 7556–7580. [Google Scholar] [CrossRef] [PubMed]
  45. Zgarbová, M.; Otyepka, M.; Šponer, J.; Hobza, P.; Jurečka, P. Large-scale compensation of errors in pairwise-additive empirical force fields: comparison of AMBER intermolecular terms with rigorous DFT-SAPT calculations. Physical Chemistry Chemical Physics 2010, 12, 10476–10493. [Google Scholar] [CrossRef] [PubMed]
  46. Robustelli, P.; Piana, S.; Shaw, D.E. Developing a molecular dynamics force field for both folded and disordered protein states. Proceedings of the National Academy of Sciences 2018, 115, E4758–E4766. [Google Scholar] [CrossRef] [PubMed]
  47. Martín-García, F.; Papaleo, E.; Gomez-Puertas, P.; Boomsma, W.; Lindorff-Larsen, K. Comparing molecular dynamics force fields in the essential subspace. PLoS One 2015, 10, e0121114. [Google Scholar] [CrossRef] [PubMed]
  48. Fotopoulos, V.; Grau-Crespo, R.; Shluger, A. Thermodynamic analysis of the interaction between metal vacancies and hydrogen in bulk Cu. Physical Chemistry Chemical Physics 2023, 25, 9168–9175. [Google Scholar] [CrossRef]
  49. Bodlos, R.; Fotopoulos, V.; Spitaler, J.; Shluger, A.; Romaner, L. Energies and structures of Cu/Nb and Cu/W interfaces from density functional theory and semi-empirical calculations. Materialia 2022, 21, 101362. [Google Scholar] [CrossRef]
  50. Grau-Crespo, R.; Hamad, S.; Catlow, C.R.A.; De Leeuw, N. Symmetry-adapted configurational modelling of fractional site occupancy in solids. Journal of Physics: Condensed Matter 2007, 19, 256201. [Google Scholar] [CrossRef]
  51. Jay, A.; Gunde, M.; Salles, N.; Poberžnik, M.; Martin-Samos, L.; Richard, N.; de Gironcoli, S.; Mousseau, N.; Hémeryck, A. Activation–Relaxation Technique: An efficient way to find minima and saddle points of potential energy surfaces. Comp. Mater. Sci. 2022, 209, 111363. [Google Scholar] [CrossRef]
  52. Woodley, S.M. Knowledge Led Master Code Search for Atomic and Electronic Structures of LaF3 Nanoclusters on Hybrid Rigid Ion–Shell Model–DFT Landscapes. The Journal of Physical Chemistry C 2013, 117, 24003–24014. [Google Scholar] [CrossRef]
  53. Larsen, A.H.; Mortensen, J.J.; Blomqvist, J.; Castelli, I.E.; Christensen, R.; Dułak, M.; Friis, J.; Groves, M.N.; Hammer, B.; Hargus, C.; others. The atomic simulation environment—a Python library for working with atoms. Journal of Physics: Condensed Matter 2017, 29, 273002. [Google Scholar]
  54. Mishin, Y.; Mehl, M.; Papaconstantopoulos, D.; Voter, A.; Kress, J. Structural stability and lattice defects in copper: Ab initio, tight-binding, and embedded-atom calculations. Physical Review B 2001, 63, 224106. [Google Scholar] [CrossRef]
  55. Zhang, B.; Hu, W.; Shu, X. Theory of embedded atom method and its application to materials science—atomic scale materials design theory. Hunan University Publication Press, Changsha, China 2003.
  56. Zhang, J.M.; Wen, Y.N.; Xu, K.W. Calculation of the formation energies of isolated vacancy and adatom–vacancy pair at low-index surfaces of fcc metals with MAEAM. Applied surface science 2007, 253, 3779–3784. [Google Scholar] [CrossRef]
  57. Jin, H.S.; An, J.D.; Jong, Y.S. EAM potentials for BCC, FCC and HCP metals with farther neighbor atoms. Applied Physics A 2015, 120, 189–197. [Google Scholar] [CrossRef]
  58. Wen, Y.N.; Zhang, J.M. Surface energy calculation of the fcc metals by using the MAEAM. Solid State Communications 2007, 144, 163–167. [Google Scholar] [CrossRef]
  59. Kühne, T.D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V.V.; Seewald, P.; Stein, F.; Laino, T.; Khaliullin, R.Z.; Schütt, O.; Schiffmann, F.; others. CP2K: An electronic structure and molecular dynamics software package-Quickstep: Efficient and accurate electronic structure calculations. The Journal of Chemical Physics 2020, 152, 194103. [Google Scholar] [CrossRef]
  60. Perdew, J.P.; Chevary, J.A.; Vosko, S.H.; Jackson, K.A.; Pederson, M.R.; Singh, D.J.; Fiolhais, C. Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation. Physical review B 1992, 46, 6671. [Google Scholar] [CrossRef]
  61. Lippert, B.G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. Molecular Physics 1997, 92, 477–488. [Google Scholar] [CrossRef]
  62. Fletcher, R. Practical methods of optimization; John Wiley & Sons, 2013.
  63. Broyden, C.G. A class of methods for solving nonlinear simultaneous equations. Mathematics of computation 1965, 19, 577–593. [Google Scholar] [CrossRef]
  64. VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. The Journal of chemical physics 2007, 127, 114105. [Google Scholar] [CrossRef]
  65. Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Physical Review B 1996, 54, 1703. [Google Scholar] [CrossRef]
  66. Wu, X.; You, Y.W.; Kong, X.S.; Chen, J.L.; Luo, G.N.; Lu, G.H.; Liu, C.; Wang, Z. First-principles determination of grain boundary strengthening in tungsten: Dependence on grain boundary structure and metallic radius of solute. Acta Materialia 2016, 120, 315–326. [Google Scholar] [CrossRef]
  67. Scheiber, D. Segregation and embrittlement of gold grain boundaries. Computational Materials Science 2021, 187, 110110. [Google Scholar] [CrossRef]
  68. Bodlos, R.; Scheiber, D.; Spitaler, J.; Romaner, L. Modification of the Cu/W Interface Cohesion by Segregation. Metals 2023, 13, 346. [Google Scholar] [CrossRef]
  69. Campbell, G.H.; Belak, J.; Moriarty, J.A. Atomic structure of the σ5 (310)/[001] symmetric tilt grain boundary in tantalum. Scripta materialia 2000, 43, 659–664. [Google Scholar] [CrossRef]
  70. Campbell, G.; Belak, J.; Moriarty, J. Atomic structure of the Σ5 (310)/[001] symmetric tilt grain boundary in molybdenum. Acta materialia 1999, 47, 3977–3985. [Google Scholar] [CrossRef]
  71. Momma, K.; Izumi, F. VESTA: a three-dimensional visualization system for electronic and structural analysis. Journal of Applied Crystallography 2008, 41, 653–658. [Google Scholar] [CrossRef]
  72. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool. Modelling and Simulation in Materials Science and Engineering 2009, 18, 015012. [Google Scholar] [CrossRef]
  73. Hirel, P. Atomsk: A tool for manipulating and converting atomic data files. Computer Physics Communications 2015, 197, 212–219. [Google Scholar] [CrossRef]
  74. Henkelman, G.; Uberuaga, B.P.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. The Journal of chemical physics 2000, 113, 9901–9904. [Google Scholar] [CrossRef]
  75. Pulay, P. Convergence acceleration of iterative sequences. The case of SCF iteration. Chemical Physics Letters 1980, 73, 393–398. [Google Scholar] [CrossRef]
  76. Farrow, M.; Chow, Y.; Woodley, S. Structure prediction of nanoclusters; a direct or a pre-screened search on the DFT energy landscape? Physical Chemistry Chemical Physics 2014, 16, 21119–21134. [Google Scholar] [CrossRef] [PubMed]
  77. Lazauskas, T.; Sokol, A.A.; Woodley, S.M. An efficient genetic algorithm for structure prediction at the nanoscale. Nanoscale 2017, 9, 3850–3864. [Google Scholar] [CrossRef] [PubMed]
  78. Barkema, G.; Mousseau, N. The activation–relaxation technique: an efficient algorithm for sampling energy landscapes. Computational materials science 2001, 20, 285–292. [Google Scholar] [CrossRef]
  79. Barkema, G.; Mousseau, N. Event-based relaxation of continuous disordered systems. Physical review letters 1996, 77, 4358. [Google Scholar] [CrossRef] [PubMed]
  80. Salles, N.; Richard, N.; Mousseau, N.; Hémeryck, A. Strain-driven diffusion process during silicon oxidation investigated by coupling density functional theory and activation relaxation technique. The Journal of Chemical Physics 2017, 147, 054701. [Google Scholar] [CrossRef] [PubMed]
  81. Ganchenkova, M.; Yagodzinskyy, Y.; Borodin, V.; Hänninen, H. Effects of hydrogen and impurities on void nucleation in copper: simulation point of view. Philosophical Magazine 2014, 94, 3522–3548. [Google Scholar] [CrossRef]
  82. Zhou, W.H.; Zhang, C.G.; Li, Y.G.; Zeng, Z. Creeping motion of self interstitial atom clusters in tungsten. Scientific reports 2014, 4, 1–4. [Google Scholar] [CrossRef]
  83. Carlberg, M.; Münger, E.; Hultman, L. Dynamics of self-interstitial structures in body-centred-cubic W studied by molecular dynamics simulation. Journal of Physics: Condensed Matter 2000, 12, 79. [Google Scholar] [CrossRef]
  84. Seeger, A.; Mehrer, H. in: Vacancies and Interstitials in Metals, Eds. A. Seeger, D. Shumacher, W. Schilling, and J. Diehl, North-Holland Publ. Co., Amsterdam 1970 (p. 1).
  85. Martínez, E.; Uberuaga, B.P. Mobility and coalescence of stacking fault tetrahedra in Cu. Scientific Reports 2015, 5, 1–5. [Google Scholar] [CrossRef]
  86. Ackland, G.J.; Vitek, V. Many-body potentials and atomic-scale relaxations in noble-metal alloys. Physical review B 1990, 41, 10324. [Google Scholar] [CrossRef]
  87. Wang, H.; Rodney, D.; Xu, D.; Yang, R.; Veyssiere, P. Pentavacancy as the key nucleus for vacancy clustering in aluminum. Physical Review B 2011, 84, 220103. [Google Scholar] [CrossRef]
  88. Osetsky, Y.N.; Barashev, A.; Zhang, Y. On the mobility of defect clusters and their effect on microstructure evolution in fcc Ni under irradiation. Materialia 2018, 4, 139–146. [Google Scholar] [CrossRef]
  89. Matsukawa, Y.; Zinkle, S.J. One-dimensional fast migration of vacancy clusters in metals. Science 2007, 318, 959–962. [Google Scholar] [CrossRef] [PubMed]
  90. Fitzgerald, S. Crowdion–solute interactions: Analytical modelling and stochastic simulation. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 2015, 352, 14–17. [Google Scholar] [CrossRef]
  91. Derlet, P.M.; Nguyen-Manh, D.; Dudarev, S. Multiscale modeling of crowdion and vacancy defects in body-centered-cubic transition metals. Physical Review B 2007, 76, 054107. [Google Scholar] [CrossRef]
  92. Hehenkamp, T.; Berger, W.; Kluin, J.E.; Lüdecke, C.; Wolff, J. Equilibrium vacancy concentrations in copper investigated with the absolute technique. Physical Review B 1992, 45, 1998. [Google Scholar] [CrossRef] [PubMed]
  93. Razumovskiy, V.; Divinski, S.; Romaner, L. Solute segregation in Cu: DFT vs. Experiment. Acta Materialia 2018, 147, 122–132. [Google Scholar] [CrossRef]
  94. Wurmshuber, M.; Burtscher, M.; Doppermann, S.; Bodlos, R.; Scheiber, D.; Romaner, L.; Kiener, D. Mechanical performance of doped W–Cu nanocomposites. Materials Science and Engineering: A 2022, 857, 144102. [Google Scholar] [CrossRef]
  95. Ebner, A.S.; Jakob, S.; Clemens, H.; Pippan, R.; Maier-Kiener, V.; He, S.; Ecker, W.; Scheiber, D.; Razumovskiy, V.I. Grain boundary segregation in Ni-base alloys: A combined atom probe tomography and first principles study. Acta Materialia 2021, 221, 117354. [Google Scholar] [CrossRef]
  96. Huang, Z.; Chen, F.; Shen, Q.; Zhang, L.; Rupert, T.J. Uncovering the influence of common nonmetallic impurities on the stability and strength of a Σ5 (310) grain boundary in Cu. Acta Materialia 2018, 148, 110–122. [Google Scholar] [CrossRef]
  97. Huang, Z.; Chen, F.; Shen, Q.; Zhang, L.; Rupert, T.J. Combined effects of nonmetallic impurities and planned metallic dopants on grain boundary energy and strength. Acta Materialia 2019, 166, 113–125. [Google Scholar] [CrossRef]
  98. Suzuki, A.; Mishin, Y. Interaction of Point Defects with Grain Boundaries in fcc Metals. Interface Science 2003, 11, 425–437. [Google Scholar] [CrossRef]
  99. Han, J.; Thomas, S.L.; Srolovitz, D.J. Grain-boundary kinetics: A unified approach. Progress in Materials Science 2018, 98, 386–476. [Google Scholar] [CrossRef]
Preprints 70620 g001
Figure 1. (a) Relative energies of the lowest energy configurations of V 3 , V 4 , V 5 , and V 6 clusters identified using the SOD code and optimized with DFT. The highest energy configuration for each cluster is used as a reference. (b) Binding and formation energies per added vacancy calculated using EAM and DFT. For each cluster, the lowest energy configuration identified through DFT is considered. (c) Schematic illustrations of the lowest energy cluster configurations for V 3 - V 6 clusters optimized using DFT. The configurations 1–5 are shown in the same energetic order as in (a). Cu atoms are blue, vacancies are red, whereas SIAs are shown in white.
Figure 1. (a) Relative energies of the lowest energy configurations of V 3 , V 4 , V 5 , and V 6 clusters identified using the SOD code and optimized with DFT. The highest energy configuration for each cluster is used as a reference. (b) Binding and formation energies per added vacancy calculated using EAM and DFT. For each cluster, the lowest energy configuration identified through DFT is considered. (c) Schematic illustrations of the lowest energy cluster configurations for V 3 - V 6 clusters optimized using DFT. The configurations 1–5 are shown in the same energetic order as in (a). Cu atoms are blue, vacancies are red, whereas SIAs are shown in white.
Preprints 70620 g001

(a) (b)
Preprints 70620 g002
Figure 2. (a) Migration barriers of V 1 , V 2 and V 3 calculated using DFT and CI-NEB. The most favorable migration path is included for each cluster. (b) Schematics of the most favorable migration paths of V 2 and (c) V 3 identified using DFT/NEB. Black arrows the direction of vacancy displacement during cluster migration. Both clusters diffuse via a single vacancy hop mechanism.
Figure 2. (a) Migration barriers of V 1 , V 2 and V 3 calculated using DFT and CI-NEB. The most favorable migration path is included for each cluster. (b) Schematics of the most favorable migration paths of V 2 and (c) V 3 identified using DFT/NEB. Black arrows the direction of vacancy displacement during cluster migration. Both clusters diffuse via a single vacancy hop mechanism.
Preprints 70620 g002
Preprints 70620 g003
Figure 3. (a) Barriers for preferred paths of V 4 and V 5 cluster migration calculated using DFT and ARTn. For both V 4 and V 5 , the initial and final configurations of the migration path correspond to the already identified lowest energy structures. (b) Local minima configurations for the migration of (i) V 4 and (ii) V 5 , with the same order as seen in (a). Both clusters migrate via a single vacancy hop mechanism.
Figure 3. (a) Barriers for preferred paths of V 4 and V 5 cluster migration calculated using DFT and ARTn. For both V 4 and V 5 , the initial and final configurations of the migration path correspond to the already identified lowest energy structures. (b) Local minima configurations for the migration of (i) V 4 and (ii) V 5 , with the same order as seen in (a). Both clusters migrate via a single vacancy hop mechanism.
Preprints 70620 g003
Preprints 70620 g004
Figure 4. (a) (i) [100](210) Σ 5 twin grain boundary used as an example of a vacancy sink. Different colors are used for atoms that sit in different planes. Numbers indicate the examined positions of a vacancy in the GB. The red-colored region illustrates the interaction range between GBs and vacancies. (ii) 304-atom (left) and 912-atom (right) Σ 5 GB simulations cells used for DFT and EAM, respectively. (b) Single point (prior to relaxation) and fully relaxed absorption energies of mono-vacancies at different sites in the grain boundaries. (c) Absorption energy calculations of fully relaxed mono-vacancies at varying distances from the GB. The interaction range between vacancies and GB (absorption region) is approximately 6.9 Å wide. In both (b) and (c) plots, inset images of the grain boundaries are included.
Figure 4. (a) (i) [100](210) Σ 5 twin grain boundary used as an example of a vacancy sink. Different colors are used for atoms that sit in different planes. Numbers indicate the examined positions of a vacancy in the GB. The red-colored region illustrates the interaction range between GBs and vacancies. (ii) 304-atom (left) and 912-atom (right) Σ 5 GB simulations cells used for DFT and EAM, respectively. (b) Single point (prior to relaxation) and fully relaxed absorption energies of mono-vacancies at different sites in the grain boundaries. (c) Absorption energy calculations of fully relaxed mono-vacancies at varying distances from the GB. The interaction range between vacancies and GB (absorption region) is approximately 6.9 Å wide. In both (b) and (c) plots, inset images of the grain boundaries are included.
Preprints 70620 g004
Preprints 70620 g005
Figure 5. (a) (i–v) Lowest energy initial configurations of di-vacancies prior to relaxation. Configurations are shown from (i) to (v) in energetic order, with (v) resulting in the lowest energy GB configuration after relaxation. Both DFT and EAM identified the configuration shown in image (v), where two vacancies relax at a 3.7 Å distance, as the most favorable initial di-vacancy configuration. (b) The lowest energy fully-relaxed configuration with a di-vacancy absorbed by the grain boundary using EAM. The di-vacancy causes the shift of the main symmetry axis of GB. (c) The lowest energy fully-relaxed configuration of the same GB with absorbed di-vacancy using DFT.
Figure 5. (a) (i–v) Lowest energy initial configurations of di-vacancies prior to relaxation. Configurations are shown from (i) to (v) in energetic order, with (v) resulting in the lowest energy GB configuration after relaxation. Both DFT and EAM identified the configuration shown in image (v), where two vacancies relax at a 3.7 Å distance, as the most favorable initial di-vacancy configuration. (b) The lowest energy fully-relaxed configuration with a di-vacancy absorbed by the grain boundary using EAM. The di-vacancy causes the shift of the main symmetry axis of GB. (c) The lowest energy fully-relaxed configuration of the same GB with absorbed di-vacancy using DFT.
Preprints 70620 g005
Preprints 70620 g006
Figure 6. (a) Energy profiles for migration of a di-vacancy from the bulk towards the [100](210) Σ 5 grain boundary axis calculated using both EAM and DFT with different cell sizes. (b) Illustrates the di-vacancy configurations along the migration path from the bulk into the GB. Configurations are shown in the same order as in (a). In the final configuration, the two vacancies are separated by 3.7 Å along the projected axis.
Figure 6. (a) Energy profiles for migration of a di-vacancy from the bulk towards the [100](210) Σ 5 grain boundary axis calculated using both EAM and DFT with different cell sizes. (b) Illustrates the di-vacancy configurations along the migration path from the bulk into the GB. Configurations are shown in the same order as in (a). In the final configuration, the two vacancies are separated by 3.7 Å along the projected axis.
Preprints 70620 g006
Table 1. Calculated migration barriers of V 1 –V 6 clusters.
Table 1. Calculated migration barriers of V 1 –V 6 clusters.
Cluster Barrier Method
V 1 0.65 DFT/NEB
V 2 0.40 DFT/NEB
V 3 0.52 DFT/NEB
V 4 0.84 DFT/EAM/ARTn
V 5 0.84 DFT/EAM/ARTn
V 6 0.96 EAM/ARTn
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
bsky
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2026 MDPI (Basel, Switzerland) unless otherwise stated