QM-accurate simulation of dislocation motion in γ-Ni and α-Fe using a hybrid multiscale approach

We present an extension of the ‘learn on the fly’ method to the study of the motion of dislocations in metallic systems, developed with the aim of producing information-efficient force models that can be systematically validated against reference QM calculations. Nye tensor analysis is used to dynamically track the quantum region centered at the core of a dislocation, thus enabling quantum mechanics/molecular mechanics simulations. The technique is used to study the motion of screw dislocations in Ni-Al systems, relevant to plastic deformation in Ni-based alloys, at a variety of temperature/strain conditions. These simulations reveal only a moderate spacing (∼5 Å) between Shockley partial dislocations, at variance with the predictions of traditional molecular dynamics (MD) simulation using interatomic potentials, which yields a much larger spacing in the high stress regime. The discrepancy can be rationalized in terms of the elastic properties of an hcp crystal, which influence the behavior of the stacking fault region between Shockley partial dislocations. The transferability of this technique to more challenging systems is addressed, focusing on the expected accuracy of such calculations. The bcc α-Fe phase is a prime example, as its magnetic properties at the open surfaces make it particularly challenging for embedding-based QM/MM techniques. Our tests reveal that high accuracy can still be obtained at the core of a dislocation, albeit at a significant computational cost for fully converged results. However, we find this cost can be reduced by using a machine learning approach to progressively reduce the rate of expensive QM calculations required during the dynamical simulations, as the size of the QM database increases.


I. INTRODUCTION
The need to produce accurate dynamical representations of chemical processes involving defects in metallic systems has been steadily increasing in recent years.An accurate description of dislocations, including plastic deformation and the effects of impurity atoms, could provide crucial and currently missing insight for industrially relevant problems, for example helping to elucidate the effect of Rhenium in Ni-based superalloys [1], or the fundamental mechanisms underlying hydrogen embrittlement of steels [2,3].Density functional theory (DFT) is a tremendously useful tool and nowadays a standard technique in many branches of physics, chemistry, and materials science [4].State-of-the-art implementations of this technique are, however, limited to a few hundred atoms (growing to a few thousand for linear scaling approaches [5], although these are currently not suitable for modeling metallic systems) and to a ∼10 ps timescale.Much larger model systems are required to accommodate the elastic field of dislocations or to describe systems closely resembling experimental samples [6].These systems are usually modelled using interatomic potentials whose accuracy has been the subject of significant research effort.
A wide range of different interatomic potentials for fcc and bcc metals have been independently developed in recent years, all sharing the idea that the strength of a bond cannot be described solely using a pair interaction model.To accurately model metallic bonding, a correction term needs to be included, depending on the local environment of the atoms and, in particular, on the site density, defined as the sum over neighboring sites of a cohesive potential in the seminal work by Finnis and Sinclair (FS) [7].While the original FS model is still very relevant for this class of materials [8], many refinements of this concept have been proposed.These include, but are not limited to, (i) the embedded atom model (EAM), originally proposed by Daw, Murray, and Baskes for the study of crystal defects in metallic systems [9] and widely applied for molecular dynamics simulations of pure or binary systems (e.g., fracture in α-Fe [10] or plastic deformation at the γ /γ interface in Ni alloys [6]), (ii) the modified EAM (MEAM), where angular dependence of the bonding is explicitly considered, allowing deviations from purely metallic bonding to be described [11], (iii) the reactive force field (ReaxFF) model, originally developed for modeling hydrocarbons but extended over the years to a number of different systems, including fcc and bcc metals [12], (iv) the charge optimized many-body (COMB) potential, capable of capturing charge transfer effects and recently employed for the study of misfit dislocation at the γ /γ interface in Ni-based superalloys [13].Unfortunately, the high accuracy of these classes of interatomic potential is generally limited to a restricted range of atomic species and phases.The development of accurate parametrized potentials including even more complex chemical environments poses significant challenges and is difficult to systematically validate, in particular for use in truly predictive situations [14].In practice, for many (but not all) "chemomechanical" problems involving realistic metallic systems the use of classical MD is not viable, as suitably general and accurate ("reactive") force fields are not available, nor is it clear how to produce fitting databases a priori guaranteed to contain the information necessary to describe all the chemical processes which might be encountered during dynamics.A range of embedding approaches have been proposed to address this length scale problem.Pioneer work by Rao, Woodward et al. [15] demonstrated how it is possible to embed a DFT calculation in a much larger elastic medium.This approach has been applied to compute strengthening effects of impurities in alloys [16] and core structures in Ni-based superalloys [17].The QM-CADD approach [18] allows dislocations to be seamlessly transported from atomistic to continuum regions.Multiprecision QM/MM techniques [19] allow chemically complex regions to be described with QM accuracy, while faraway atoms provide the correct elastic embedding, typically accurately described by MM modeling.This is particularly useful for the modeling of chemomechanical processes, in which a chemical process in a small region affects the macroscopic response of a system.Within the framework of energy-based QM/MM techniques, the total energy of a system is written as the sum of a classical term for the MM region, a quantum term for the QM one and an interaction term.This artificial interface might lead to an incorrect charge density and to spurious ionic forces permeating into the QM region.This issue has been recently addressed by constraining the charge density at the QM/MM interface to an accurate target density that reflects the atomistic configuration within the MM region [20].The method reproduces the QM charge density and the magnetic moments of bulk materials (Al, Fe, Pd) and produces a reasonable edge dislocation core structure for bulk iron.In another work, the QM region is embedded in a suitable atomistic virtual environment (not the MM region), so that the QM calculation is performed in a periodic system rather than a cluster of atoms, thus reducing the effect of Friedel oscillations propagating from the free surface [21].These QM/MM approaches provide accurate structural models that can be used for identifying the equilibrium configuration of crystal defects in metals with QM accuracy.However, none of these approaches can model temperature or free-energy effects.
The "learn on the fly" (LOTF) scheme [22] is a QM/MM algorithm based on the fitting of a corrective energy function, fitted to reproduce target QM forces in the core region and augmented by a predictor/corrector algorithm for computational efficiency.This method has been successfully adopted to simulate a number of problems related to fracture in semiconductor materials, including impurity-driven scattering mechanisms [23] and stress corrosion [24].Recently, it has been shown that the LOTF approach can be applied to the study of dislocations in metallic systems [25] and can be extended to compute potential energy barriers to dislocation glide [26].Moreover, dynamical QM/MM simulations produce large amounts of first-principles data, which can be used as part of a training database for a machine learning (ML) model aimed at the construction of an accurate force field (FF), capable of replacing the QM method whenever "new" configurations (not well represented within the training database) are not encountered.Recent years have seen a growing research interest in such data-driven approaches with the goal of significantly increasing the time scales accessible to simulations.There are two broad classes of ML approaches: (i) learn "once and for all" approaches where the potential energy surface of the target systems is fitted once-and-for-all using, e.g., neural networks [27] or Gaussian process regression [28], with the latter recently applied to iron [29][30][31]; (ii) "on-the-fly learning" approaches where the force model is trained adaptively using a database of QM forces calculated during an MD run [32][33][34].The second class of method has also lead to the development of an ML extension of the LOTF algorithm [32].
The present work is structured as follows.In Sec.II an overview of the theoretical methods is provided, focusing on the main differences between previous implementations of the LOTF algorithm, aimed at the simulation of fracture in brittle materials [23], and the present one, capable of tracking a dislocation core under high temperature conditions.These simulations, reported in Sec.III and performed over a range of temperature and strain conditions, and for increasing chemical complexity (i.e., increasing Al percentage within the γ matrix), reveal an improved description of dislocation core geometry over traditional MM simulations.In particular, the distance between Shockley partial (SP) dislocations is shown to have a weak dependence on the temperature/strain conditions of the simulation.A different behavior is observed for classical MM simulations, in particular under high stress conditions.This discrepancy can be rationalized as a poor EAM description of elasticity in hcp Ni/Al crystals.In Sec.IV we report future prospects for the application of the LOTF scheme in metallic systems.First, we address the feasibility of extending the use of this QM/MM scheme to α-Fe.We find that the method remains sound and sufficient precision can be attained, but the necessary QM calculations become significantly more challenging because of the complex electronic structure of iron, connected with its magnetic properties.Finally, we show that state-of-the-art ML methods can in principle achieve accurate prediction of QM forces for the materials and systems presented in this work and could thus provide appropriate support to the use of LOTF techniques.

A. Simulation cell
The Ni/Al system is modelled using a quasi-2D periodic cell comprised of a quadrupole of screw dislocations in the γ matrix, as illustrated in Fig. 1.The system is oriented with z orthogonal to the dislocation line ξ = 1 √ 2 [ 110] and y normal to the close packed family of planes {111}.The x axis, parallel to [11 2], is the direction of motion of a gliding a 2 110 {111} dislocation.The system length along the directions orthogonal to ξ is ∼170 Å, corresponding to a distance of ∼80 Å between screw dislocations.This maximum distance between dislocation cores is on average preserved during a MD simulation due to the symmetry of the system.The system length along z is chosen to be 5 Å (two elementary units), periodic along the dislocation line, effectively decreasing the number of atoms required for the QM treatment and thus increasing the overall simulation speed.We used the linear elastic displacement field solution to add four perfect screw dislocations to the matrix and then relaxed them with the EAM potential, leading to dissociation into pairs of Shockley partials.We selected this geometry principally to enable direct comparisons between the QM/MM and EAM models; we note that the choice of a quadrupole dislocation arrangement does not fully decouple the dislocations from one another [35,36].Before carrying out MD simulations, the perfect crystal is initially rescaled according to the expected thermal expansion for the EAM potential, as calculated in Ref. [37].

B. The QM/MM scheme
First-principles calculations are performed using the projected-augmented wave (PAW) implementation of the Vienna ab initio simulation package (VASP) [38,39].The Perdew, Burke, and Ernzerhof (PBE) functional [40] is used for the exchange-correlation term.Brillouin zone integration is performed with Monkhorst-Pack grids [41].The electronic smearing is included following the Methfessel-Paxton scheme [42], with a broadening of 0.1 eV.We follow the LOTF simulation scheme, as described in Ref. [23], modifying the buffer termination and the QM selection scheme, as detailed below.The saturation of broken bonds with hydrogen atoms, which previously proved to be essential to obtain accurate forces on QM clusters in covalent materials [23], is not necessary for metals, due to the lack of directionality of bonds.Accurate DFT forces for Ni clusters can instead be obtained with a bare surface providing a sufficiently large buffer region is used, as discussed in our earlier work [25].The QM selection scheme mirrors the topology and structure of the material defects.The topological criteria used for studying fracture in brittle materials are not appropriate here, as a clear difference in the coordination number at the dislocation core is not always observed as the position of the dislocation evolves at finite temperatures.Instead, we used the Nye tensor α to identify and track dislocation cores.This was originally developed in Ref. [43] and more recently extended to atomistic systems.It is now well established as a tool capable of identifying dislocation cores by analyzing the local environment of an atom [44] and it is here used to locate the QM region.We note that the average force error of the EAM potential with respect to DFT is proportional to the magnitude of the Nye tensor components, as demonstrated in Ref. [25].In the specific case of a screw dislocation dissociated into SPs, only two components of this tensor are relevant, α 33 and α 31 , with the former contributing to the screw component of the Burgers vector and the latter to the edge.The atomic-resolved screw component α 33 is equivalent for the two SPs and integrating it over the plane normal to the dislocation line produces the total Burgers vector.In the case of the edge component α 31 , contributions from the two SPs have opposite signs, and the surface integral vanishes.However, integrating only in the neighborhood of the SP core would produce the modulus of the edge component of the SP dislocation.A plot of these quantities if shown in Fig. 10.
At first sight, the edge component of the Nye tensor seems more suitable than the screw one to locate an SP during a molecular dynamics simulation, as the difference in sign would always make the two SPs distinguishable.Some complications, however, arise when a system at finite temperature is considered, primarily caused by the short distance between SPs in γ -Ni.While the distance between SP centers is 5 Å, each core has a finite size of about 2 Å, corresponding to a localized region (∼1 Å) between the SPs in which the Nye tensor is vanishing (or at least negligible when compared to the cores).The introduction of temperature in the system causes these areas to overlap from time to time, resulting in destructive interference in the edge components of the Nye tensor, which vanishes by definition for the case of perfect screw dislocation.This problem is not encountered for the screw component, as the interference is in this case constructive, and the two QM regions would simply collapse into a single one until the next dissociation event, which would be accounted for by recalculation of the Nye tensor, as implemented in our LOTF algorithm.
A key advantage of buffered QM/MM schemes such as the LOTF is the possibility of updating the selection of the QM atoms during a simulation, thus allowing a moving defect to be tracked using a QM cluster of minimal size [19].In this work, we initially identify the QM region with the position of two atoms giving large contributions to the Nye tensor (α 33 > 0.04 Å −1 ), with a separation d comparable with the equilibrium SP separation value of 5 Å.Two circular regions of radius d/2 are selected, and the position of the cores is updated as the average of the atomic positions using α 33 as weight.Note that this quantity does not need to be updated at every time step of the molecular dynamics simulation (2 fs), as the thermal oscillation frequency of the core is smaller than that of individual atoms.A QM calculation is similarly not necessary at every step, and the force error can be controlled using a predictor-corrector scheme, as implemented in the LOTF method [22].Testing reveals that a QM calculation is necessary every 7 time steps (14 fs) to maintain the force error below 0.1 eV/Å.The QM region is updated using a hysteretic selection algorithm [19], with inner and outer selection radii of d in = 3 Å and d out = 5 Å, respectively, to reduce fluctuations in the set of QM atoms.These small QM regions are enough to ensure good accuracy, as atoms farther than 4 Å from dislocation cores exhibit bulklike behavior [25].
The buffer atoms are required to decouple the QM region from the surface states, thus ensuring accurate QM forces for the QM/MM embedding.In the case of γ -Ni, this is obtained with a 5 Å buffer [25], leading to clusters containing ∼50 atoms.Due to the geometry of the simulation cell, these clusters are periodic along the line direction z.Vacuum is introduced only in the (x,y) plane and the Brillouin zone sampled with a 1 × 1 × 7 MP grid which reduces to 4 k points once time reversal symmetry is taken into account.
The force convergence threshold of 0.1 eV/Å is chosen as a compromise between the accuracy of the simulation and its computational cost.This value is compatible with the typical EAM force errors for bulklike systems with respect to a DFT reference, as shown in Ref. [25].More accurate results can be achieved with the LOTF method, albeit at larger computational cost.For example, errors below 0.05 eV/Å on an individual QM atom can be achieved using a cluster of ∼250 atoms, corresponding to an increase in computational cost by a factor of 5 3 = 125 for cubic scaling DFT (see also Appendix C).The full system is shown in Fig. 1, coloring the atoms according to the screw Nye tensor component α 33 , and with the QM selection algorithm illustrated in the right panel.

C. Machine learning approach
The redundancy in the configurations visited during the course of an ab initio MD simulation allows the construction of an accurate FF directly built from the QM data.Gaussian process (GP) regression provides a parameter-free data interpolation technique [45,46] that can be efficiently employed to make use of these data.A reference database is built by associating the forces {f i } N i=1 on a collection of N atoms with the corresponding local environments {ρ i } N i=1 of the same atoms.Such local environments are typically stored as the positions (relative to the central atom) of all the particles within an appropriate cutoff radius.The measure of the correlation between forces associated with different configurations ρ, ρ is defined to be a matrix-valued kernel function K(ρ, ρ ).This can also be augmented by encoding specific prior knowledge, e.g., the symmetries of the system [30].The GP prediction of the force associated with a new configuration ρ * is given by evaluating its correlation with all the elements in the dataset [45,46], where the vectors {c i } can be found in closed form [47,48].
There are various reasons why GP regression can be preferable to standard fitting of a fixed function of interatomic positions.For instance, since no system-specific functional form has to be identified, this kind of regression generally requires very little ad hoc tuning-model (kernel) selection can be easily streamlined.Furthermore, GPs typically improve their predictive accuracy as the database size grows, so that an arbitrary accuracy can be achieved when enough data is available for a suitable choice of kernel functions.
To test the capabilities of ML models, we collected a database by performing ab initio MD simulations of nickel and iron in their crystalline cubic phases, using 4 × 4 × 4 supercells of their respective fcc and bcc elementary cells.The MD simulations were carried out within the canonical ensemble using a Langevin thermostat with a friction coefficient γ = 2 ps −1 .Brillouin zone sampling was carried out using a 4 × 4 × 4 MP grid, and symmetries in both real and reciprocal space were disregarded to improve the simulation accuracy; other DFT parameters are as described in Sec.II B above.This database was then randomly separated into a training set and a test set.To build a kernel function incorporating the symmetries of the system we first constructed a permutationally invariant distance as described in Refs.[30,31,49].We then used this distance in a squared exponential kernel [45], symmetrized over the point group of the crystal (for full details refer to Ref. [30]).

III. LOTF SIMULATION OF DISLOCATIONS IN γ NICKEL
We investigated the glide of dislocations in γ -Ni at a range of temperatures, stresses, and alloy compositions, using the previously detailed implementation of the LOTF scheme.The main difference expected to be found between the MM and QM/MM simulations is the local ordering of the atoms in the proximity of the dislocation core.Discrepancies of this kind may lead to incorrect MM predictions for the separation distance between SPs, as observed in former multiprecision works, see, e.g., Ref. [50].We define the equilibrium separation distances between SP cores as the average of the atomic positions within half the stacking fault length of each core weighted by their respective Nye tensor screw components.Other approaches that have been previously used include measuring the length of the hcp stacking fault region using common neighbor analysis (CNA) [51]; this is always larger as it includes contributions from the two core regions.We prefer the Nye-based definition as it remains stable at high temperatures and also has the advantage for QM/MM purposes of centering the core on the atoms for which the deviation between DFT and EAM forces is maximal (as shown in Ref. [25]).For the EAM-relaxed configuration, we obtain an SP core separation of 5 Å with the Nye method and 10 Å with the CNA method, consistent with previous literature results using the same potential in Ref. [51].We note that this low temperature value for the separation is likely to be correct within the EAM scheme, as it only depends on the shear modulus and on the ISF energy, which are both accurately reproduced by the potential [37].The validity of this description as a function of temperature, stress, and chemical complexity is addressed in this work.Note that applying a shear strain to the system is not expected to affect the SP separation distance, with any deviation from this behavior representing an inadequacy of the modeling of interatomic interactions in the system.The applied shear stress also determines the velocity of a gliding dislocation.We have decided to use finite stress values in all our simulations, as this would allow for extensive testing of the QM region tracking algorithm presented in this work.
We initially thermalized the system at a target simulation temperature for 5 ps, using the EAM potential for all atoms.After thermalization, a shear strain yz is applied, and the simulation is carried on for another 1 ps.A further 0.5 ps of QM/MM dynamics follows, to re-equilibrate the system after the change of the Hamiltonian, before production calculations.Production QM/MM calculations were carried out for 2 ps.The MM and the QM/MM calculations are found to predict different core geometries and separations for the SP dislocations.The SP separation distances, averaged over time and over the four dislocations in the simulation cell, are reported in Table I for systems under different conditions of temperature, stress, and chemical composition.In the low stress, low temperature (50 K, 100 MPa) regime the MM and the QM/MM results are completely consistent, with no significant difference observed in the 0% to 5% Al concentration interval.This is also consistent with the accuracy of EAM-predicted stacking-fault energies and elastic moduli.A good agreement between the two methods is also found in the high temperature (1000 K) regime, as long as the percentage of Al in the matrix is small (5% or below).For a larger percentage of Al (15%), the distance between SPs increases significantly, consistently with the overestimation of the ISF energy incurred by the EAM compared with the DFT predictions (see Appendix A).In these conditions, we observe the EAM-modelled pair of partial dislocations develops a very large separation distance of ∼11 Å, which fluctuates but never reverts to the initial value in the simulation time scale.The corresponding QM/MM simulation does not exhibit deviations from the equilibrium value of ∼5 Å, as reported in Fig. 2. The corresponding initial and final configurations are also shown in the figure, coloring atoms by Nye tensor to highlight the position of the dislocation cores.We next consider the high-stress low-temperature case (300 MPa, 50 K).In the low Al composition (5%) case, the distances measured with QM/MM are still broadly in line with the previously observed results, whereas the MM ones are larger by a factor of two.
In the yet more complex case of high Al composition (15%) under high stress, the average MM-predicted distance between SPs is 15 Å, the largest average value encountered in this work.We interpret this to be due to the combined effects of a "soft" hcp phase and underestimated ISF for high Al compositions.Under these conditions, the QM/MM separation is larger than that observed in the previous cases.It is thus likely that under these working conditions larger QM regions, covering the whole stacking fault, have to be used in order to avoid spurious SP geometries.
Inspection of the trajectories reveals that large SP separation distances are in most cases caused by a failure in the description of the leading/trailing dislocation pair.The trailing dislocation, in particular, fails to follow the leading one until the latter has moved several steps, thus resulting in unphysically large separation distances.We also note that in the case of high temperature simulations, the MM separation distance is often shorter than the QM/MM one, in spite of its average value being larger.This often leads to constriction and cross-slip processes, not observed under the QM/MM descriptor.Selected trajectory frames are reported in Fig. 3, showing an example case.Dislocation constriction is observed after 0.6 ps in the MM simulation.At 0.8 ps the dislocation, still constricted, starts to cross slip.A perfect screw dislocation glide is then observed, until another cross-slip process occurs at 1.4 ps.The dislocation dissociates then into SPs in a plane equivalent to the starting one.At this point, the distance between partials becomes larger, also due to the presence of a larger number of Al atoms in the current layer.In the QM/MM case, the distance between SP oscillates, but cross-slip phenomena are never observed.
To summarize, we find that the EAM potential gives reasonable results only at low temperature, for small values of the applied deformation and only when low concentrations of impurities are present in the matrix.As previously noted, all the inconsistencies between QM and MM simulations derive from an overly soft description of the HCP phase, making it easier to stretch (large stress) or contract (high temperature) depending on the simulation parameters.Accounting for this issue could greatly improve the accuracy of the specific EAM potential under examination.With this in mind, we dedicate the remainder of this section to the identification of the discrepancy between the QM and the QM/MM approaches.This can be separated into two independent contributions arising from the elastic energy of the hcp phase and from the chemical composition of the alloy.The latter is not difficult to rationalize: The Al-Ni pair interaction is fitted to the γ phase.This description is not transferable to the Al impurity in the Ni matrix, due to the similarity between the chemical environment of these atoms, and the binding energy of the defect is overestimated, leading to a smaller ISF energy.A more complete discussion of this issue is provided in Appendix A. The strain component of the EAM error is more subtle.We started by noticing that, while the energy difference between the fcc and the hcp phases is included in the fitting database for the EAM potential, the elastic constants for the latter are not.We calculated and compared the elastic energies using an orthorhombic simulation cell, with the z axis parallel to the [111] direction.The hcp and fcc crystals unit cells contain four and six atoms, respectively.These were periodically reproduced along z to obtain two easily comparable, stoichiometric structures containing six {111} layers each.Compressive and tensile strains in the range ±10% were applied along z with 1% increments, to compute energy differences with respect to the equilibrium reference, displayed in Fig. 4. Good agreement between DFT and EAM is observed for fcc, consistent with the potential having been fitted to the elastic constants [37].On the other hand, the elastic energy of the hcp crystal is underestimated.This is particularly true in the near anharmonic regime for compressive strains, with errors larger than 25% for strains above 5%.This last feature is surprising, as expectations of the limited transferability of a classical potential would normally suggest a different likely scenario, i.e., one with similar predicted elastic energetics in two materials for which the first shell of neighbors of each atom are the same.We therefore attempted to trace the origin of this effect, considering the interatomic potential in more detail.Within the EAM model, the total energy of a system of N atoms is represented as the sum of a pair interaction term and an embedding many-body one: where r i j are the interatomic distances.The pair interaction φ(r) is a generalization of the Lennard-Jones potential.The embedding energy term F (ρ) is constructed using a density function ρ(r), taken to be a combination of exponentials, and obtained by fitting the EAM to an equation of state [37].
The elastic energy associated with the embedding term is illustrated in Fig. 4(b) and is practically identical for the two structures.This is fully consistent with the exponential decay of the contributions to the electronic density energy term, which is the ansatz of the EAM approach.The deviation between the two crystals, leading as discussed above to an underestimation of the elastic energy in the hcp structure, must therefore originate in the pair interaction term, as confirmed by Fig. 4(a).This can be further understood by considering the contribution to the pair energy term of individual shells of neighbors.For clarity, the hcp and fcc structures are illustrated in Fig. 5, indicating atoms from different shells with different colors.The reference central atom is shown in black, and atoms in the same (111) plane are neglected, as they are equivalent in the two crystals.The first two shells of neighbors (in yellow and orange) contain six atoms at 2.49 Å and six atoms at 3.52 Å in both the fcc and hcp structures.The equilibrium distances of the third and fourth neighbor shells are different in the two structures.The distance values are provided in Fig. 5 and reported in Fig. 6(a) as color-filled circles (green and blue for fcc, red and teal for hcp) placed on top of the pair interaction plot (dashed black line).The pair interaction contribution for the elastic energy of these second two shells is further explored in Figs.6(b) and 6(c).In the case of tension, the major contribution for the fcc structure (6 × 6.9 meV) is associated with the six atoms of the third neighbor shell [green in Fig. 6(c)].The hcp structure contains only two atoms in the third neighbor shell of the reference atom, located in second-neighbor (111) planes above and below.These atoms [red in Fig. 6(b)] produce a scaled down total contribution of 2 × 6.7 meV.The difference from the fcc case is not balanced by the contribution of the more populated 12-atom hcp fourth shell [teal in Fig. 6 contributes 12 × 1.2 meV.In the case of compression, the (teal) fourth shell is associated with a significant contribution of 12 × −3.4 meV, which is this time more than enough to balance the effect of the (green) third neighbor fcc shell (6 × −6.7 meV).A further stabilizing contribution to the energy in the hcp crystal is then associated with the presence of the (red) third neighbor two-atom shell (2 × −3.4 meV).The fourth neighbor fcc shell [blue in Fig. 6(c)] does not meanwhile provide a significant stabilizing contribution, again because these atoms lay close to the pair interaction cutoff distance.As a result of adding these negative contributions to the energy upon compression, a soft hcp phase is predicted, which accounts for the discrepancy of the curves in Fig. 5(b).This observation illustrates the challenges of capturing the elastic response of both the fcc and hcp by means of a single pair potential form and provides a practical example of the subtle transferability issues that interatomic potentials may encounter, even state of the art high-quality ones that very accurately and usefully describe the phase they were originally designed for, as is the case here.

A. Extension to other materials
The investigations above exemplify how possible shortcomings of interatomic potentials are not generally easy to identify a priori, as is clear from Sec. III, while in many cases they can be a posteriori rationalized by comparison with QM-based methods.With this in mind, we find that the QM/MM embedding method provides a useful tool for validating potentials, as well as an alternative route to carry out MD simulations when interatomic potentials are not sufficiently accurate.The remaining challenge in the latter scenario is that QM/MM-based MD simulations, even when time accelerated by using a predictor/corrector recipe as in the LOTF scheme, are much slower than classical MM calculations.This is especially significant for materials where the supporting QM calculations require large cluster convergence radii due to the basic physics of the target system (or to its chemical complexity, e.g., due to alloying).In these cases, it becomes mandatory to minimize the number of necessary QM calculations.To exemplify this situation and illustrate a possible way forward (using a data-driven approach-see the next subsection), we describe preparatory work for the simulation of a screw dislocation in α-Fe in a multiprecision scheme.This complements earlier results for nickel already presented in Ref. [25].
Iron is a challenging material for QM/MM approaches, as the magnetization is much larger than in nickel alloys.Spurious spin rearrangement at the surfaces of the QM clusters will be correspondingly higher, so that Fermi-level crossing effects become more likely to propagate to the center and affect the calculated forces, making large QM clusters necessary.This is explicitly shown in Appendix C, by plotting the projected density of states of a single QM atom for spherical clusters of increasing size, comparing the cases of iron and nickel.Various heuristic approaches to mitigate the slow convergence of the central atom to a bulklike electronic structure with respect to the width of the buffer are investigated here, including (i) imposing a larger Fermi-energy smearing width and (ii) allowing the QM cluster to host a net charge or (iii) constraining the magnetization at the surface to its equilibrium value.The usage of dipole corrections is found to not affect the computed forces significantly, indicating that electric dipoles generated on the QM clusters by surface effects are negligible.We studied the convergence of the QM forces with respect to cluster radius for two systems: a 4 × 4 × 4 supercell of bulk α-Fe (128 atoms) and a larger one centered on a a 2 111 {110} screw dislocation.Classical calculations were performed using a high quality EAM potential [52].The simulated systems were thermalized at 1200 K by means of classical MD simulations, following the same protocol detailed in Ref. [25].For each of the 10 independent snapshots considered, clusters of increasing size were carved by considering all the neighbors within a certain distance from the test QM atom.The convergence of the QM forces in ferromagnetic bulk iron is displayed in Fig. 7(a).These results indicate that forces fall within our chosen 0.1 eV/Å accuracy threshold with clusters containing ∼110 atoms or more, while the corresponding nickel system reached the same accuracy level with only 55 atoms.For an interesting comparison, when auxiliary calculations on the same Fe clusters in fictitiously imposed non-spin-polarized conditions are carried out [teal symbols in Fig. 7(a)], the convergence of the QM forces is much faster.It should be noted that the forces obtained from these nonmagnetic DFT calculations deviate significantly from those obtained in the ferromagnetic case.These results clearly indicate that force convergence is coupled to the magnetization of the Fe physical systems.To further investigate this effect, we calculate the atom-resolved magnetic moment of 0 K α-Fe clusters, containing full neighbor shells with a well-defined number of atoms, using Bader charge analysis [53,54] on DFT charge densities.The obtained values are shown in Fig. 7(b) as a function of the distance from the test atom, for clusters of different sizes.The deviation of the magnetic moment for the Fe atom at the center of a QM cluster from its bulk value of 2.15 μ B can be as high as 20% in a cluster containing three neighbor shells ["cl3" orange symbol in Fig. 7(b)].Moving to atoms neighboring the surface, the magnetization can be as much as 50% larger than its bulk value.We find that, for the central atom and its nearest neighbors to be unaffected by surface effects, large clusters including as many as twelve neighbor shells must be used.
Next, we tried constraining the magnetic moment of the atoms in the QM clusters to its bulk value to mimic the average environment far from surfaces, in order to investigate whether this could yield converged QM forces for smaller cluster sizes.This was done by computing the magnetic moment from a Mulliken analysis in a sphere centered on each atom and adding a penalty term to the total energy proportional to the square of the deviation of each magnetic moment from its target bulk value, summed over all atoms in the cluster.The strength of the constraint can be tuned by varying the penalty term prefactor λ, which ideally should be as large as possible to force a nearly constant magnetization across the QM cluster but still small enough not to induce instabilities in the self-consistent electronic structure optimization loop.Our atomic force results for clusters extracted from the 1200 K trajectory and containing ∼60 atoms [the same used for the calculations reported in Fig. 7(a)] are presented in Table II.These calculations show larger deviations of the forces from TABLE II.Relation between the strength of the magnetic moment constraint and the force error.
f (eV/Å) 0.09 ± 0.02 0.10 ± 2 0 .12± 0.03 the correct QM target values as the λ parameter is increased from zero (unconstrained calculation) to higher values, indicating that constraining magnetic moments cannot easily be used to speed up force convergence.We next investigated the convergence of QM forces with respect to the size of the buffer for atoms neighboring a dislocation core in α-Fe.The resulting data are reported in Fig. 7(c).As for screw dislocations in γ -Ni, producing reference DFT calculations for the whole (∼50 000 atom) system is not feasible, so the results of DFT calculations performed on a large (∼200 atom) cluster are taken as the reference.Here our results suggest that forces on atoms at the core of a screw dislocation are on average converged within 0.1 eV/Å if the clusters used in a QM/MM scheme contain ∼140 atoms.We find that this conclusion is robust to changes in the Fermi energy smearing width.Large smearing widths produce a more even balance in the occupancy of Fermi level states, such as the spurious cluster surface states which might conceivably hinder the force convergence on central atoms if unevenly occupied.However, repeating calculations using a 0.5 eV smearing width does not yield significantly improved convergence.This can be appreciated in Fig. 7(c), where these results are presented together with the original curve obtained with a 0.1 eV smearing width and express deviation from the same reference force (200 atoms, 0.1 eV smearing width).Comparing the force error plots obtained for different values of the smearing width suggests that the intrinsic error introduced by using a higher smearing may balance any beneficial smoothing effects of the kind discussed above for system sizes at which the absolute error is still above the target threshold.Since larger clusters are still required to reduce the average force error to values lower than 0.1 eV/Å, using large smearing widths brings no convergence benefit.Similarly, adding or removing electrons in a self-consistent fashion does not yield any significant improvement in the convergence of the QM forces, so that achieving occupancy balance of quasidegenerate surface cluster states by charge addition or depletion (in the ±5 e range, while using a neutralizing background charge density in the DFT calculations) does not yield improved force convergence with QM cluster size.This is shown in Fig. 7(d), for a 59 atom cluster at 1200 K, showing an essentially flat force error profile.
Taken together, these results suggest that force convergence to the accuracy required to make QM/MM calculations of dislocation motion practical can be attained in unconstrained local spin density calculations on α-Fe systems of interest.While this is good news, it is clear that very large QM-zone clusters are required to fully decouple the central atom from surface states in QM/MM calculations aiming at a quantitative convergence to their reference (while practically unattainable) fully-QM results.To enable practical simulations we next note that the overall cost of the calculations could be still very significantly reduced by judicious storage and reuse of the information produced.Namely, machine learning (ML) techniques could then be used to predict atomic forces using a database of reference atomic configurations and their associated accurate QM forces.Ideally, for any new atomic configuration "similar enough" to previously encountered ones, a ML interpolation scheme would be sufficient to predict the forces on atoms necessary to continue a MD simulation, so that only genuinely new configurations would ever require new QM calculations [32].The accuracy of ML methods based on Gaussian process (GP) regression are tested in the next subsection.

B. Gaussian process regression for bulk metals
We tested our GP models (see Sec. II C) on sample trajectories extracted from ab initio MD simulations of nickel and iron in their crystalline cubic phases.Figure 8 shows scatter plots of DFT force components versus those predicted by both EAM [37,55] and the GP models.These plots reveal that EAM predictions (blue points) are adequate on average in the region of configuration space where the forces are small.As the magnitude of the forces increases, the EAM predictions incur significant systematic errors with respect to the DFT reference forces (overestimated in Ni and slightly underestimated in Fe: this feature depends however on the exact choice of the potential).The accuracy of the GP predictions (red points), on the other hand, does not depend on the modulus of the force being predicted, an effect that can be related to the absence of a fixed functional form systematically "biasing" the predictions.
To illustrate how the ML force field can be systematically improved, we generated a second ab initio dataset for iron, selecting atoms neighboring an isolated vacancy.The mean GP error on forces as a function of training points in the database is shown in Fig. 9.This reveals that the average EAM accuracy can be improved by using moderately sized training sets (∼100 configurations) and that the target 0.1 eV/Å force-error threshold is also easily reachable for the systems considered.We note here that the computational cost of any new GP prediction increases with the number of training points [cf.Eq. ( 1)].However, since the scaling is only linear, the GP model remains orders of magnitude faster than computing forces by direct DFT calculations even for the largest database used in this work (N = 320 points).These results (as well as those of many recent works on the construction of ML force fields for a variety of materials [30,31,[56][57][58]) suggest that the production of data-driven force fields could allow for a drastic reduction of computational effort while maintaining QM accuracy.However, when simulating nonequilibrium chemomechanical phenomena, training a ML force field "once and for all" before the simulation is carried out will often not represent a sufficiently good strategy, since the system is likely to explore regions of configuration space not well represented in the initial training database.An "on-the-fly learning" scheme might be of help in these circumstances in as much as the training database could be modified or augmented adaptively.This approach represents an active field of research [32,34,59], and further developments are surely necessary to fully exploit its capabilities.Nevertheless, the good performance of GP-based force models in the present class of systems and the increasing availability of QM databases provide a promising practical background for information efficient on-the-fly learning techniques.

V. CONCLUSION
We have presented an extension to the "learn on the fly" method, aimed at the QM/MM modeling of dislocation motion in metallic bulk systems.The moving QM region (dislocation core) is tracked using an algorithm based on the Nye tensor analysis and automatically updated during the course of the QM/MM simulation, thus allowing for the selection of a (mobile) QM region of minimal size, reducing the computational cost.This approach is then used for studying the glide of a screw dislocation in γ -Ni, focusing on the local arrangement of atoms at its core and, in particular, on the length of the stacking faulted region separating Shockley partial dislocations.The study is repeated under different simulation conditions including temperature, applied shear strain, and chemical composition.The EAM-based interatomic potential, used for the MM part of this work, is found to be in good agreement with the QM/MM model at the low stress/low chemical complexity regime, even at high temperatures.Simulations under more extreme conditions reveal the development of a large separation distance between SPs during a MM simulation, behavior that is not observed using the QM/MM method.This deviation is rationalized in terms of the lack of transferability of the EAM potential, incapable of correctly describing an Al impurity in γ -Ni or the elastic response of hcp nickel.While this established the need of validation and the usefulness for this purpose of reference QM embedding methods, it comes at a significant computational cost that increases further with the complexity of the system, as seen in the case of α-Fe.However, machine learning methods appear fit to learn QM forces accurately in metallic systems for the class of problems addressed here, suggesting that the efficiency issue can be addressed using databases of atomic configurations for which QM-accurate forces are known.In particular, tests on bulk systems of γ -Ni and α-Fe reveal that DFT data can be reproduced using databases of modest size.High accuracy is obtained also in the vicinity of a point defect (vacancy).These results indicate that a QM/MM implementation is a viable tool to investigate dislocation motion in metallic systems, making it possible to accurately describe the nonequilibrium atomic configurations found in the neighborhood of dislocation cores.This scheme could be supported by a ML algorithm capable of reproducing QM-accurate forces in all encountered configurations similar to those contained in the available database, for which interpolation is all that is required, limiting the need for developing new QM information to genuinely new configurations encountered during the system's dynamical evolution.This opens the path to applying an on-the-fly learning (ML-based LOTF [32]) QM/MM schemes in these systems, allowing bridging between the typical simulation times accessible to ab initio and classical molecular dynamics.
where E tot is the total energy of the relaxed defected system and E ref the energy of a Ni crystal of the same size.The chemical potential μ(X) is evaluated from the per-atom energy of fcc crystals.The EAM value for the formation energy deviates from the DFT reference, as reported in Table III.This discrepancy can be explained by noting that the EAM potential is fitted to the energetics of the γ phase, in which the first shell of neighbors is the same as for the Al impurity.
As the interactions beyond nearest neighbor are small for the EAM potential, this constrains the energy prediction for the impurity.This can readily be seen by comparing the cohesive energy of a substitutional Al atom in γ -Ni with that of an Al atom in γ , using both the full EAM potential and a reduced version for which only nearest neighbor interactions are considered.This gives almost identical results (Table III).
The lower substitutional formation energy discussed above is a feature depending mostly on nearest neighbor interactions, and it is therefore expected to apply also to the hcp case, leading once again to the prediction of a larger distance between SPs (note that in this case this is true also for the equilibrium 0 K geometry).This is verified here by calculating the ISF energies using DFT calculations and comparing them to EAM results.The crystal is oriented along [  operation disrupts the periodicity along z, 15 Å of vacuum are introduced in this direction.This stacking configuration is then relaxed using the FIRE algorithm with a maximum force tolerance of 0.05 eV/Å, constraining atomic positions to change only along z.The ISF energy can at this point be obtained as the difference between the final configuration and the starting (fcc) one, normalized over the area of a {111} layer.We repeated this calculation for systems containing increasing amounts of Al (0%, 12.5%, 25%) in the faulted {111} layer, each time placing the Al atoms so that they were never, before or after the shift, nearest neighbors, as this would be energetically unfavorable.Our results, reported in Table IV, confirm the expected trend.As this suggests the equilibrium SP separation will also be overestimated in the presence of Al impurities, consistent with the behavior of dislocation core pairs observed in MD simulations in Sec.II B. temperature (1200 K), and we investigate the convergence of electronic properties for increasing widths of the buffer region.The same test is applied to both γ -Ni and α-Fe.The projected densities of states (PDoS) for these systems are calculated at the point using the Methfessel-Paxton smearing.For some of the smaller systems, we have verified that a more dense sampling of the BZ would not yield to significantly different results.A bulk reference is provided, using a 4 × 4 × 4 reciprocal space grid and the same smearing type and parameter adopted for the clusters.
The results for a selection of clusters of different size are displayed in Fig. 11 with colored lines; the bulk reference is indicated with a gray filled curve.The PDoS of the cluster systems are more localized with respect to the bulk reference, and a very sharp peak can be observed for the smallest cluster examined.A large number of atoms (≈300 and ≈400 atoms for Ni and Fe, respectively) are required to properly converge the electronic structure.Starting from clusters of 176 atoms, the PDoS for the Ni atom becomes similar to the bulk reference, with the exception of the major peak at 0.5 eV below the Fermi level, where the height is overestimated.Clusters of at least 300 atoms are required to mitigate this effect.A different deviation from the bulk reference is observed for these system: The bulk reference is indicated with a gray filled curve.The configuration is extracted from a MD trajectory at 1200 K.The Fermi level is fixed to 0 and indicated with a vertical black solid line.The positive and negative portions of the vertical axis denote spin-up and spin-down contributions to the pDoS, respectively.a minor intermediate peak between the major one at 0.5 eV and Fermi.
In the case of iron, the convergence of the electronic properties is more challenging, in particular due to the several crossings of the Fermi level observed for the spin-up channel [see for instance the 260 atoms cluster in Fig. 11(a)].The major peak of bulk Fe, at −0.75 eV, is split into a series of minor peaks, and clusters as large as 330 atoms still present this feature.Starting from clusters of 414 atoms the PDoS starts to resemble the reference profile.Note that each of the reported PDoS exhibits a different value at Fermi.These large variations are not observed in nickel, for which the Fermi level is close to extrema of the PDoS (a minimum for the spin-up channel and a minimum for the spin down) and exhibits therefore small fluctuations with respect to the accuracy of the calculation.Another way to investigate the convergence of electronic properties is monitoring how the charge density at the central atom and its space integral evaluated by means of Bader analysis converge with cluster size.As charge rearrangements are not expected in this system, we have opted for studying the convergence of the magnetic moment.We present these results in Fig. 12(b), together with the convergence of the QM force on the same atom, shown in panel (a).Similar oscillations of the magnetic moment can be observed for both γ -Ni and α-Fe.Larger errors are observed in the case of iron, consistent with the larger intrinsic magnetization of this material.Notably, most cluster calculations for iron would predict a magnetic moment larger than the bulk one.This is consistent with the very large magnetization observed for the surface atoms (see main text).The 0.1 eV/Å force accuracy threshold is reached with the first (≈40 atoms) and the fourth (≈140) atom shell for Ni and Fe, respectively.This makes a QM/MM treatment of nickel much cheaper than for iron, as discussed in the main text.However, a force accuracy as low as 0.05 eV/Å is eventually reached for both systems using QM clusters larger than 250 atoms.Note however that force convergence does not correspond to a full convergence of all electronic properties, as is, e.g., clear from the DoS plots in Fig. 11 from the oscillating magnetic moment errors reported in Fig. 12(b) [see also panel (b) in Fig. 7, main text].These results indicate that the QM clusters used for practical LOTF calculations do not describe the electronic properties of the system with sufficient accuracy.This justifies the usage of "tricks" for accelerating the convergence of QM forces, including large values of Gaussian width for electronic smearing approaches.

FIG. 1 .
FIG. 1. Quadrupole of screw dislocations in γ -Ni during a QM/MM simulation at room temperature.Atoms are colored according to the value of the screw component of the Nye tensor, which was also used to define the QM region.For clarity, only one of the two overlapping QM regions is shown in the lower right panel.

FIG. 2 .
FIG. 2. Time-averaged SPs separation and geometries for a QM/MM and a MM simulation at high temperature conditions (1000 K, 200 MPa) for an alloy containing 15% Al.Al atoms are represented using larger spheres.

FIG. 3 .
FIG. 3. Comparison between molecular dynamics trajectories inthe MM and QM/MM approach.Selected frames from simulation at (1000 K, 200 MPa, and 5% Al).Al atoms are represented using larger spheres.

FIG. 4 .
FIG.4.Elastic energy due to strain along the [111] direction calculated using DFT (blue lines) and EAM (red lines) for nickel fcc (a) and hcp (b) crystals.Associated contributions to the EAM elastic energy: pair interaction (c) and embedding energy (d) for nickel fcc (red circles) and hcp (blue squares) crystals.

FIG. 6 .
FIG.6.Pair interaction potential for Ni atoms.Panel (a): full energy profile (dashed line) and neighbors positions (circles).The number of atoms in each shell is also indicated.In panels (b) and (c) the pair interaction is plotted for the third and the fourth shells of neighbors, which are different in the two structures.The equilibrium (full circles) and deformed (5% strain, empty circles) positions are indicated, together with the elastic energy variation per atom upon deformation, for each shell of neighboring atoms.

FIG. 7 .
FIG. 7. (a) Convergence of DFT forces with respect to the size of the cluster for ferromagnetic (FM) and artificial nonmagnetic (NM) bulk iron structures.(b) Magnetization error with respect to bulk for atomic shells in differently-sized iron clusters.(c) Convergence of the DFT forces on atoms neighboring a dislocation core.(d) Effect of adding/removing charge from the clusters.

FIG. 8 .
FIG. 8. Scatter plots showing accuracy of GP regression vs EAM on 200 randomly selected testing configurations.A 160-configurations training database was used for the GP models.

FIG. 9 .
FIG. 9. Average force error for Fe atoms within 3 Å of a vacancy as a function of the size of the database.

FIG. 10 .
FIG. 10.Screw (a) and edge (b) components of the Burgers vector for the relaxed MM core structure of screw dislocations dissociated into Shockley partials.

FIG. 11 .
FIG. 11.Projected density of state of the central atom (QM region) of a spherical buffer of increasing width for α-Fe (a) and γ -Ni (b).The bulk reference is indicated with a gray filled curve.The configuration is extracted from a MD trajectory at 1200 K.The Fermi level is fixed to 0 and indicated with a vertical black solid line.The positive and negative portions of the vertical axis denote spin-up and spin-down contributions to the pDoS, respectively.

FIG. 12 .
FIG.12.Force error (a) and magnetic moment error (b) with respect to the bulk reference for a monoatomic QM region for clusters of increasing size in γ -Ni (blue) and α-Fe (red).

TABLE I .
Average SP separation at finite temperature for different stress regimes and compositions, as evaluated by Nye tensor analysis of MM and QM/MM simulations.
The molecular dynamics trajectories created during this research are openly available from the University of Warwick Research Portal at [60].

TABLE III .
Substitutional formation energy of an Al impurity in γ -Ni evaluated using EAM and DFT (a).Cohesive energies and for an Al atom in the γ phase, as calculated with the full EAM potential and with a truncated version limited to nearest neighbors interactions.The latter produces by construction the same results in the (undefected) γ phase (b).For instance, they can lower the ISF energy, thus increasing the equilibrium SPs separation.Rhenium impurities are, e.g., very useful for the design of efficient Ni superalloys, due to their beneficial effects on tertiary creep deformation (the so-called rhenium effect).Here we investigate the energetics of the Al substitutional impurity in γ , for which we find a less satisfactory agreement with DFT predictions for the EAM implementation.The model systems used in DFT calculations are constructed by substituting a Ni atom with an Al in a 3 × 3 × 3 cubic cell containing 108 bulk atoms.Details of the calculations are the same as reported in Sec.II B above.Geometries are optimized within a 1 meV total energy tolerance, using a conjugate gradients algorithm and 4 × 4 × 4 MP k-point mesh for BZ sampling.The substitutional formation energy is evaluated as sub 1 00], [11 2], and [111].A 2 × 2 × 6 supercell is used (18 {111} layers).The stacking fault is generated by shifting the nine upper planes by b/ √ 3, where b is the nearest neighbor distance.Since this

TABLE IV .
Intrinsic stacking fault energies at different Al concentrations.