Modelling defects in Ni–Al with EAM and DFT calculations

We initially test the EAM potential against DFT data, to verify its reliability in describing the properties of Nickel alloys. We use an EAM model developed by Mishin in 2004 [23], generally considered to be amongst the most reliable interatomic potentials for modelling $\gamma /{{\gamma}^{\prime}}$ interfaces and plasticity-related phenomena in Ni-based alloys. This is due to its accuracy in reproducing the correct lattice parameters for both phases as well as quantities directly related to dislocations such as the intrinsic stacking fault energy of Nickel. A new EAM potential for simulating Ni–Al alloys was made available in 2009 [24]. It is however less accurate for describing $\gamma /{{\gamma}^{\prime}}$ interface geometries due to its underestimation of the lattice parameter for ${{\gamma}^{\prime}}$. The potential is fitted to the properties of the B2 NiAl phase, and the only quantity related to ${{\gamma}^{\prime}}$ included in the potential is its formation energy.

To gain more insight in the issue of transferability of EAM potentials in Ni systems, we preliminarily compare DFT and EAM results for a selection of structural properties. For a more complete study and a comparison with experimental data, refer to [23]. Our electronic structure calculations have been performed using the VASP package [25], adopting the PBE/GGA approximation for the exchange and correlation energy [26]. Brillouin zone (BZ) sampling was carried out using Monkhorst–Pack (MP) grids [27] with different meshes for different simulation cells, as specified for each case below. We used a 300 eV plane-wave energy cutoff, the projector augmented wave (PAW) method [28], and Fermi energy smearing with a Gaussian (Methfessel–Paxton) broadening of 0.1 eV [29].

As a first step, we evaluated the equilibrium lattice parameter for both the γ and ${{\gamma}^{\prime}}$ phases using four-atom cubic elementary cells for different volumes (11 configurations, ±5% from experimental value). A $13\times 13\times 13$ MP grid is used for BZ sampling. The resulting energies were fitted to the Murnaghan equation of state. Both DFT and EAM predict equilibrium lattice parameters of 3.52 Å and 3.57 Å for γ and ${{\gamma}^{\prime}}$ respectively. Bulk moduli obtained from fitting are reported in table 1.

Table 1. Elastic properties (in GPa) for bulk models of Ni and $\text{N}{{\text{i}}_{3}}\text{Al}$ systems. Comparison between DFT and EAM data.

	DFT (γ)	EAM (γ)	DFT $\left({{\gamma}^{\prime}}\right.$)	EAM (${{\gamma}^{\prime}}$)
B	192	181	177	176
C11	268	237	236	236
C12	161	151	147	154
C44	127	125	126	127

Elastic moduli are important quantities for a QM/MM simulation, as any mismatch between the two descriptors would result in a discontinuity at the interface between the classical and the quantum regions [17]. The ability to reproduce the correct elastic constants (once the lattice constants are aligned by parameter tuning) is thus an important feature for an interatomic potential. We evaluated the elastic constants with a linear fit to the stress/strain curve, using stress tensors calculated using the virial (derivative of potential energy with respect to cell vectors). For these calculations, we use a 400 eV plane-wave cutoff energy, necessary for obtaining accurate diagonal elements of the stress tensor for the ${{\gamma}^{\prime}}$ phase. The resulting data are presented in table 1, and show good agreement between the DFT and EAM approaches.

Vacancies are crucial defects to enable plasticity behaviour through dislocation climb. The EAM potential used here was explicitly fitted to reproduce the experimental values for vacancy formation and diffusion energy [23]. It is therefore important to check the consistency of these values with the DFT predictions. We used a $3\times 3\times 3$ supercell of the conventional 128 atom cubic system for the γ phase, with one atom removed to model an isolated vacancy. The system was relaxed to its structural minimum within a 1 meV energy tolerance using a conjugate gradient algorithm and $4\times 4\times 4$ MP $\mathbf{k}$-point grid. We then calculated the vacancy formation energy and migration barrier, the latter using the NEB method, which involved the optimisation of 3 intermediate images by a RMM-DIIS quasi-Newton, variable metric algorithm for the QM calculations and the FIRE algorithm [30] for the MM calculations. Our calculated results for the γ phase are reported in table 2. The agreement between EAM and DFT is overall good, with the interatomic potential overestimating the reference DFT values properties by only around 10%.

Table 2. Formation and migration energies for a vacancy and intrinsic stacking fault energy per unit surface in the γ phase of Ni.

	DFT	EAM
$E_{f}^{v}\left(\text{eV}\right)$	1.45	1.57
$E_{m}^{v}\left(\text{eV}\right)$	1.08	1.19
${{E}_{\text{ISF}}}\left(\text{eV}~\text{n}{{\text{m}}^{-2}}\right)$	0.85	0.84

Close agreement between the calculated values of intrinsic stacking fault energy (ISF) using DFT and EAM is also important, as perfect dislocations in fcc metals are always dissociated into Shockley Partials (SP) and the distance between them is inversely proportional to the ISF energy. To investigate this issue, we constructed a model system by orienting the crystal with the z direction perpendicular to the normal to the $\left\{1\,1\,1\right\}$ close packed family of planes. The x and y directions were chosen perpendicular to the $\left[\bar{1}\,1\,0\right]$ and the $\left[1\,1\,\bar{2}\right]$ crystal directions, respectively. A twelve $\left(1\,1\,1\right)$ layer supercell was initially created as a bulk reference. We then removed one layer, and rescaled the simulation cell in order to produce a periodic model system containing a single stacking fault defect, where the distance between planes with hcp stacking and their nearest periodically repeated images is  ∼20 Å. For the structural optimisation we used a conjugate gradient algorithm. A MP grid of $19\times 11\times 3$ $\mathbf{k}$-points was used to sample the BZ for QM calculations. The ISF energy value was obtained as a direct subtraction of the bulk reference energy (renormalised for the correct number of atoms) from the energy of the defected system fully relaxed up to 1 meV tolerance, after dividing by the area of a single $\left(1\,1\,1\right)$ layer. Our results are reported in table 2, and once more indicate close match between EAM and DFT results.

The generally good agreement between EAM and DFT results for systems in their structural minima provides no indication of the performance of the EAM potential when the system is displaced from equilibrium. Perhaps the simplest test configuration for non-zero forces geometries is obtained by displacing a single atom from its equilibrium structure, and checking the consistency of the forces predicted by the two models. We therefore performed a set of calculations using a $3\times 3\times 3$ supercell with one of the atoms displaced along the $\left[1\,1\,0\right]$ direction i.e. towards one of its nearest neighbours. The forces acting on both atoms are presented in figure 1(a) as a function of the applied displacement, revealing a significant discrepancy between the two approaches when bond lengths are modified by 10% or more. A similar behaviour has been reported in [31] where the EAM ansatz of an exponentially decreasing charge density was studied for different deformations of a fcc crystal by measuring the QM charge density in a (vacant) octahedral site. This work found that when a second nearest neighbour of the vacant site moves towards it, the site loses charge in favour of the region comprised between the approaching atom and its nearest neighbours, which violates the EAM approach ansatz.

Zoom In Zoom Out Reset image size

The plot of the charge density distribution for a 12% bond deformation is displayed in figure 1(b). The density associated with isolated atoms was subtracted in this $\left(1\,0\,0\right)$ plane plot, so that the charge depletion from atomic sites and its redistribution into the interatomic space, typical of the metallic bonding, is immediately apparent. A charge transfer pattern similar to the one described in [31] is observed, indicating that the failure of the EAM to describe the charge rearrangement can lead to the system visiting unphysical configurations during a constant temperature MD simulations, or at the very least that the distribution of visited configurations may not correctly sample the canonical ensemble consistent with the reference QM energetics. While these discrepancies still occur as relatively rare events for bulk systems simulated at moderate temperatures, we find that large force errors are incurred when using the EAM model in the defected systems containing vacancies or dislocations, as described in section 3.

We next test the ability of the EAM potential to reproduce the DFT forces in configurations produced by finite temperature MD simulations, as would be necessary to correctly predict thermodynamic observables [32]. We start from simple bulk structures, and investigate the EAM force error temperature dependence. A $3\times 3\times 3$, 108 atom supercell was used for modelling both the γ and ${{\gamma}^{\prime}}$ phases. Temperature effects were introduced in the system by means of classical MD simulations carried out at fixed volumes in each case consistent with thermal expansion. The system was initially equilibrated for 5 ps, using a 2 fs time step. This was followed by a further 10 ps period during which configuration snapshots were collected at 1 ps intervals. The DFT forces were then calculated on the collected structures (using calculation parameters as in the previous sections) and compared with the classical values.

The error distribution plots produced by these calculations are presented in figure 2, and the associated average force errors reported in table 3. These results suggest that for perfect bulk systems the EAM interatomic potential is remarkably accurate, with more than 70% of the absolute errors contained below the 0.1 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ threshold for room temperature γ phase simulations. However, the average error increases for higher temperatures as presented in figures 2(a) and (b) and approximately doubles for T  =  1200 K (table 3). The error distribution of the relative errors is shown in figures 2(c) and (d), showing no appreciable difference between the plotted curves for different temperature values. This indicates that the errors incurred by the EAM predicted forces are approximately proportional to the forces, which implies an acceptable performance in simulating bulk systems even at relatively high (while still below melting) temperatures.

Zoom In Zoom Out Reset image size

Table 3. Average force component error of the EAM interatomic potential computed as deviation from the reference DFT calculations.

T (K)	300	600	900	1200
Force error γ (eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$)	$0.07\pm 0.02$	$0.10\pm 0.02$	$0.12\pm 0.03$	$0.13\pm 0.04$
Force error ${{\gamma}^{\prime}}$ (eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$)	$0.10\pm 0.02$	$0.14\pm 0.04$	$0.16\pm 0.04$	$0.18\pm 0.05$

Note: Averages over ten independent γ (Ni) and ${{\gamma}^{\prime}}$ ($\text{N}{{\text{i}}_{3}}\text{Al}$) 108 atom model system configurations, at different temperatures.

Our computed results for the ${{\gamma}^{\prime}}$ phase are similar, but reveal that here EAM is overall less accurate: only about 40% of the forces are accurate within 25% of the reference DFT values. This is probably due to a greater difficulty in describing the interaction between different chemical species, consistent with the fact that fitting accurate classical potential forms so that accuracy is preserved while the number of chemical species increases is notoriously challenging [33].

The main long-term objective of this work is to identify the working simulation parameters needed to run multiscale QM/MM simulations of plasticity-related processes in Ni-based superalloys. In these simulations there is no requirement for an EAM potential to correctly reproduce forces on atoms in all conditions. Providing that local configurations for which these potentials are accurate can be identified a priori, all that is needed is to make sure that the MM zone of the simulation remains within this 'comfort zone'. A QM description can then be used whenever more chemically complex configurations (e.g. involving bond breaking and reforming) are encountered, where the potential cannot be assumed to be accurate enough.

Accordingly, we next set out to explore the QM region requirement associated with our EAM potential in the proximity of a vacant lattice site of the γ matrix. We note that vacancies are not just the most common intrinsic point defects found in metal lattices: in the present context, at phase boundaries, they are also a necessary ingredient for the dislocation climb-assisted glide mechanisms generally thought to enable creep deformation in Ni-based superalloy blades, particularly in the low-stress, high-temperature conditions where plasticity remains mostly confined in the γ channels of the alloy microstructure. Our vacancy system was obtained from a $3\times 3\times 3$ γ phase bulk cell (108 atoms), upon removal of one atom. This system is small enough that fully converged DFT reference results are readily obtainable, and large enough to make the interaction between the vacancy and its repeated images negligible for the present purpose. As in section 2.2, temperature effects were introduced by carrying out constant temperature EAM MD simulations, this time setting the temperature to 1200 K, in line with with typical operating temperatures of superalloy turbines.

Our goal is to investigate the accuracy of the EAM potential by comparing the atomic forces it produces with those obtained from DFT calculations. If discrepancies indicate that using a QM region is required, we wish to determine both the size of this QM region and the width of the associated 'buffer' region necessary for decoupling core quantum atoms from the surface states to ensure accurate QM forces are obtained [17].

Plots of the EAM error, both absolute and relative, with respect to DFT references are reported in figures 3(a) and (c). The different line colours correspond to different distances from the vacancy: 2.5 ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ on average for nearest neighbour (n.n.) atoms and 3.5 Å for next-nearest neighbours (n.n.n). Results for atoms from the same frames located more than 4.3 Å from the vacancy, and thus expected to be 'bulk-like', are also included for reference, and for these atoms about 50% of the data points have an absolute force error below 0.1 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ and only 20% exceed 0.2 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$. Here, the error is rarely larger than 0.4 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$, in line with the results of section 2.2. For the atoms neighbouring the vacancy the force error distribution is broader. While in about 40% of the data points the error is still below 0.1 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ a much larger population of forces deviating by 0.2 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ appears (about 35%), with the maximum errors being as large as 0.8 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$. Inspecting the relative error distribution function also makes clear that the relative error for these atoms is likely to be much larger than for bulk-like atoms. The error plot measured for next-nearest neighbours show an intermediate behaviour.

Zoom In Zoom Out Reset image size

We note that the EAM potential examined here is in fact by design accurate enough for modelling the dynamics of vacancies in the γ phase of Nickel alloys, as both the formation and the migration energies are included in the fitting database used to develop the potential [23]. However, although the EAM forces integrate the correct migration barrier, we find here that the local description of the potential energy surface (of which the calculated forces are simply tangent vectors) is not as accurate for the n.n. of vacancies as it is for the bulk-like atoms. As a result, the potential could easily fail to describe the interaction between a vacancy and other defects such as dislocations, for which the potential has not been fitted, so that a more accurate treatment is required. We therefore proceed to define a quantum region surrounding vacancies for accurate QM/MM simulations.

Implementing an accurate QM/MM treatment of this system requires a separate analysis of the convergence of the quantum forces with respect to the size of the QM cluster. This is presented in figures 3(b) and (d) including results for the isolated vacancy and the bulk system for comparison. For the bulk, we started by selecting our target QM atoms to be those with forces exceeding 0.5 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ and performed DFT calculations for Ni clusters of increasing size (denoted in the figure by buffer width radius) centred on each of these atoms. The cluster radii were chosen to lie between neighbour shells to reduce number fluctuations associated with thermal disorder. The amount of vacuum surrounding clusters was set to 10 Å in all calculations, after verifying that this is sufficiently large for the interactions between clusters and their periodic images to be negligible. Force errors were averaged over ten configurations.

For systems including a vacancy calculations were performed on larger clusters centred on the defect itself. In this case, we defined our QM regions as the first shell of neighbours of the vacant sites. For each of the ten system snapshots considered, we examined QM atoms with forces exceeding 0.5 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$, and averaged our convergence results to produce the plots reported in figure 3. The QM convergence behaviour is similar for the the vacancy and bulk cases. Good convergence is already obtained when the width of the buffer zone surrounding the QM atoms in the clusters is about 4.2 Å or bigger. These values correspond to including the first three shells of neighbours of a fcc crystal in the quantum total calculation. The radius of the corresponding cluster is about 7 Å for the system centred on the vacancy. While the error is flat when further increasing the cluster radius in the bulk-like system, the last point for the vacancy system still deviates slightly from the converged value (while being below our accuracy threshold), possibly due to residual image effects in the reference periodic calculation.

As a next step, we considered a screw dislocation in the γ phase. The crystal was oriented with orthorhombic cell axes $x=\left[1\,1\,\bar{2}\right]$, $y=\left[1\,1\,1\right]$, $z=\left[\bar{1}\,1\,0\right]$, leading to a six-atom unit cell. A $48\times 34\times 8$ supercell was then formed, and an $\frac{a}{2}\langle 1\,1\,0\rangle \left\{1\,1\,1\right\}$ screw dislocation was introduced by imposing the correct elastic displacement [34]. The resulting system is periodic along the dislocation line direction (z, also parallel to the Burgers vector). The dislocation core was tracked using the Nye tensor analysis [35]. Replication of the system along the line direction z allows carrying out averages of Nye tensor components during a molecular dynamics simulation at high temperature. This is useful to average out the local thermal fluctuations, making it easier to locate the partial dislocation cores and define the (approximately) cylindrical regions where a QM treatment is necessary. The z-averaged surface integral of the tensor can also be more easily checked to add up to the relevant Burgers vector.

A vacuum region was introduced along the x and y directions to prevent dislocation emission at the stacking fault, due to imposing periodic boundary conditions Atoms within 80 Å of initial dislocation core positions were free to move, while those outside this region were kept fixed during both the relaxation and the dynamics, in order to prevent the dislocation reaching the surface. The system was relaxed using a conjugate gradient algorithm using a 0.03 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ force convergence threshold. The atomic positions were then rescaled to account for thermal expansion, and the dynamics of the system simulated in the NVT ensemble. The system was initially equilibrated for 5 ps, after which five different system snapshots were considered, taken at intervals of 5 ps to minimise correlation. Atoms were classified into three categories according to the value of the zz component of the Nye tensor:

• (i)  
  n.n. of the dislocation core: defined as the atoms for which the screw component of the Nye tensor is above 0.04 ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$. Few of these atoms are present in each $1\,1\,0$ layer, but they contribute about 70% of the integrated Burgers vector. These atoms are usually grouped together around each partial core.
• (ii)  
  n.n.n. of the dislocation core: atoms with a smaller, but non-zero (0.02 ${{\overset{\circ}{{\text{A}}}\,}^{-2}}$) screw component of the Nye tensor, typically neighbouring the atoms of the previously defined region.
• (iii)  
  Bulk atoms: all other atoms located within 40 Å from the dislocation core(s).

Force accuracy was investigated in convergence tests only for the first two categories, as the faraway atoms' main role is to provide the elastic embedding for the region surrounding the dislocation, which was assumed in all cases achievable by classical potential (EAM) modelling. Full DFT calculations cannot be performed on the whole system for the system sizes required here, to provide a reference target values. However, the convergence of QM forces near the dislocation cores can still be monitored as before, by performing calculations on clusters of increasing size. As above, the forces acting on the central atom of each cluster are the quantity of interest in each calculation, while the remaining atoms in the QM cluster just act as a buffer [17].

Our force convergence results with respect to the cluster radius are displayed in figures 4(b) and (d). The reference forces are taken from a larger cluster, not shown in the plot, including one more atomic shell than those shown (i.e. nine neighbour shells, or 177 atoms). As for the vacancy systems, we performed ten sets of calculations for atoms which are n.n. to the core and averaged over these results. The resulting average error is monotonically decreasing for increasing sizes of the system, and falls below out accuracy threshold of 0.1 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ with a cluster radius of 5.5 Å. We note here that an additional neighbour shell is required for QM/MM modelling of dislocations, compared with the bulk and vacancy cases. However, further testing shows that the force convergence can be further improved by increasing the Fermi energy smearing width from 0.1 to 0.5 eV in the DFT calculations. This reduces the force errors for smaller clusters, so that accurate QM forces can be obtained with a three shell cluster for the as for the bulk and vacancy cases.

Zoom In Zoom Out Reset image size

The results described allow to determine which atoms require QM corrections, once more by direct comparison with the EAM predictions. We picked a set of 30 atoms from each system snapshot obtained above, and performed a DFT calculation for each of them using a well converged 5.5 Å cluster radius. The deviation of the EAM forces from the DFT forces evaluated this way were analysed separately for two crystallographic directions, one corresponding to the $\left[1\,1\,\bar{2}\right]$ dislocation glide direction and the other to the $\left[1\,1\,1\right]$ climb direction. We find that the EAM potential performs reasonably well in both cases, with 40% and 35% of the forces deviating by more than 0.2 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ for the glide and climb directions, respectively (the corresponding value for the bulk is 25%). Errors are, however, higher along the direction of motion for a gliding dislocation, and the presence of vacancies possibly yielding the formation of jogs during the diffusion will very likely increase the force error along the $\left[1\,1\,1\right]$ climb direction. Moreover, in both directions, there is a significant probability of incurring force deviations as large as 0.6 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ for both the nearest neighbours and the next nearest neighbours of the dislocation core. The next nearest neighbours are associated with larger EAM force errors in the $\left[1\,1\,\bar{2}\right]$ direction, i.e. the one linking the two partial dislocations across the stacking fault region where the atomic linkage less closely reproduces the perfect crystal bonding geometry.

Modern Ni-based superalloys include different atomic species known to improve the plasticity properties [2]. In particular, the so-called 'rhenium effect' has not been fully explored and understood at the atomic scale, mainly due to the lack of accurate tools for describing both the chemical bonding and the elastic properties for systems larger than a few hundred atoms. One of the biggest limitations of EAM (and generally) potentials is their lack of straightforward transferability to chemical environments not fully represented in the fitting database, for instance when describing a different phase or when introducing a new atomic species in the simulated system is required.

Since re-fitting a potential in such a way that its domain of application is generalised while the previous performance is in no way affected can require generating sizeable new fitting databases and is generally a non-trivial process, there is scope to inspect the convergence of QM/MM methods which can, in general, be immediately applied to new problems, due to the intrinsic generality/transferability of first-principles methods.

We thus devote this section to quantifying the size of a quantum region surrounding an impurity atom (specifically, rhenium or tungsten) in the γ phase, to determine the QM cluster radii required to obtain converged DFT forces. As for the vacancy analysis, we model the γ matrix with a $3\times 3\times 3$ supercell, this time substituting one atom with an impurity atom. Atomic positions are in all cases rescaled to match the lattice parameter predicted for pure Nickel by thermal expansion, assuming a target temperature of 1200 K. In contrast with the system preparation procedures described above where classical MD was used, here the relevant systems configurations were generated by DFT dynamical simulations, since no suitable EAM parameters are available for impurities.

The DFT calculation parameters are the same as described previously for cells of the same shape and size. The MD simulations were carried out in the NVT ensemble using a 2 fs time step and a Nosé Hoover thermostat with mass corresponding to an 80 fs oscillation period. The system was first equilibrated at 1200 K for 2 ps, and a snapshot configuration was afterwards sampled every 200 fs during the next 2 ps. In order to monitor the chemical effects affecting the QM region size region necessary in the presence of impurity atoms, we created an auxiliary system for each frame, in which the central atom was substituted back with a Nickel atom. This allows a comparison between the total charge density and magnetisation fields between the two systems.

We used Bader analysis [36, 37] to determine the atom-resolved charge and magnetisation. Our results are presented in figures 5(a) and (b). The behaviour is similar for the Re and W impurity species: the first shell of neighbours gains a small amount of charge (the net charge loss of the Rhenium atom is about 0.05 e for rhenium and 0.1 e for tungsten, shared between 12 neighbouring atoms). The second nearest neighbours lose charge to compensate for this effect, and atoms further away are unchanged (i.e. any local charge fluctuations are compatible with the temperature-induced symmetry breaking). The magnetisation of the impurity atom is zero in both cases, and the magnetic moment for the neighbouring Nickel atoms is also considerably reduced. The local magnetisation is practically converged to the bulk value (0.63 ${{\mu}_{\text{B}}}$) by the second and third shells of neighbours for rhenium and tungsten impurities, respectively. This is consistent with the larger charge transfer effects observed in the tungsten system.

Zoom In Zoom Out Reset image size

To further probe the local effects of the impurity on the surrounding metal matrix atoms, we measure the average force difference between the defected structures and the corresponding auxiliary ones. When this difference vanishes, Nickel atoms are no longer affected by the presence of the close-by defect. Our results are presented in figures 5(c) and (d). The force differences are above the 0.1 eV ${{\overset{\circ}{{\text{A}}}\,}^{-1}}$ threshold only for the first neighbouring shell, and vanish at 6 Å from the defect. These first principles results indicate that the QM zone must contain at least the first shell of neighbours for Rhenium, and possibly also the second one in model systems containing Tungsten atoms.

The standard test for determining the minimal size of the QM region involves a direct comparison between the DFT forces and those predicted by the interatomic potential, allowing to terminate the QM region as soon as the MM description is sufficiently accurate. Unfortunately, as mentioned above, no parametrisation has been made available extending this EAM potential to atomic species other than Nickel and Aluminium. An attempt to produce a ternary (Ni–Al–Re) interatomic potential based on the MEAM method for modelling processes at the $\gamma /{{\gamma}^{\prime}}$ interface has been reported in [38]. However, this ternary potential is not yet fully reliable as it has been reported to overestimate the lattice parameters and to predict spurious Rhenium clustering in the γ phase not observed by extended x-ray absorption fine structure and atom probe tomography [7, 8] and later ruled out using DFT calculations [10]. Another ternary potential was constructed by Du et al [39] based on a previous model for generic fcc structures [40]. In the present work we have tested this potential against DFT forces and found large deviations when modelling substitutional rhenium impurities in γ. This is probably due to the non-transferability of the potential (otherwise accurate in describing properties of non-defected ordered phases) to the isolated impurity system in high temperature conditions of interest here.

Once the size of the QM region containing the atomic forces that will be used in a QM/MM dynamical simulation is determined, the next step is to determine the thickness of the surrounding buffer region required to effectively decouple it from spurious surface effects. The two combined requirements will give the size (radius) of the overall cluster used in QM computations. We carried out convergence tests using clusters of different radii centred on the impurity, and the results are presented in figure 6. These reveal that forces on impurity atoms accurately converge to their reference periodic DFT values for a radius values larger than  ∼6.5 Å. Converging the QM forces on all QM region atoms requires a  ∼7.5 Å radius (∼180 atoms) cluster, while a 6.8 Å (∼130 atoms) cluster calculations already yields forces accurate enough to match the overall very good accuracy achieved by the interatomic potential in the surrounding (MM) region.

Zoom In Zoom Out Reset image size