Screw dislocation structure and mobility in body centered cubic Fe predicted by a Gaussian Approximation Potential

Abstract

The plastic flow behavior of bcc transition metals up to moderate temperatures is dominated by the thermally activated glide of screw dislocations, which in turn is determined by the atomic-scale screw dislocation core structure and the associated kink-pair nucleation mechanism for glide. Modeling complex plasticity phenomena requires the simulation of many atoms and interacting dislocations and defects. These sizes are beyond the scope of first-principles methods and thus require empirical interatomic potentials. Especially for the technological important case of bcc Fe, existing empirical interatomic potentials yield spurious behavior. Here, the structure and motion of the screw dislocations in Fe are studied using a new Gaussian Approximation Potential (GAP) for bcc Fe, which has been shown to reproduce the potential energy surface predicted by density-functional theory (DFT) and many associated properties. The Fe GAP predicts a compact, non-degenerate core structure, a single-hump Peierls potential, and glide on {110}, consistent with DFT results. The thermally activated motion at finite temperatures occurs by the expected kink-pair nucleation and propagation mechanism. The stress-dependent enthalpy barrier for screw motion, computed using the nudged-elastic-band method, follows closely a form predicted by standard theories with a zero-stress barrier of ~1 eV, close to the experimental value of 0.84 eV, and a Peierls stress of ~2 GPa consistent with DFT predictions of the Peierls potential.

Introduction

Iron and iron-based alloys are among the most important engineering materials for structural applications. Due to their widespread use under a wide range of conditions, there is a strong need to understand the atomic structures, motions, and interactions among defects, such as vacancies, interstitials, dislocations, surfaces, and interfaces.^1,2,3,4 Accurate understanding can be provided for some isolated defects using first-principles methods such as density functional theory (DFT). However, important macroscopic mechanical properties, such as yield stress, work-hardening, ductility, fracture toughness, and response to radiation exposure depend on defects’ interactions and motion. Studying defect interactions generally requires simulations at scales far larger than what is accessible using DFT with current high performance computing. Accurate interatomic potentials that enable large-scale atomistic simulations (molecular statics or dynamics) of complex defect–defect interactions under stress and at finite temperature would thus be extremely valuable. Insight from such atomistic simulations could then further guide the formulation of quantitative mechanistic theories that can predict macroscopic behavior and that are based on the relevant nanoscale defect–defect interaction mechanisms.

The mechanical strength and ductility of bcc iron are controlled primarily by the stress- and temperature-dependent mobility of screw dislocations.⁵ The mechanism of motion underlying this mobility is the nucleation of kink pairs along the straight screw dislocation, followed by lateral glide of the kinks, leading to the advance of the dislocation to the next energy minimum (Peierls valley) and the creation of a plastic slip equal to one Burgers vector.⁶ Knowing the mechanism, theories exist that can use DFT-computed inputs.⁷ However, as mentioned, direct DFT simulations of the kink-pair nucleation process are not feasible due to the sizes (tens of thousands of atoms) required for accurate modeling. Unfortunately, existing analytical interatomic potentials are unable to properly capture this crucial mechanism for the single-defect problem, nor can they reproduce underlying DFT-computable features such as the core structure and Peierls potential.^{8,9,10,11,12,13} If these features are explicitly fitted in such potentials, vibrational properties are not reproduced correctly.⁷ There are empirical potentials that show the emergence of the kink-pair nucleation mechanism, but at the same time they predict an incorrect core structure, Peierls potential and slip plane.^8,13,14,15

Quantum-mechanically based bond order potentials (BOP) have been pursued as an alternative. BOP potentials show many attractive features for Fe screw dislocations¹⁶ such as the non-degenerate core structure, the single-hump Peierls potential, and slip on {110}.¹⁷ Other features associated with screw dislocations and motion, such as the generalized stacking fault γ surfaces and non-Schmid effects were also investigated.¹⁷ However, the kink-pair nucleation process has not yet been investigated for bcc Fe with BOP potentials, although it has been studied with BOP for non-magnetic bcc systems.^18,19,20 While successful in many respects, BOP remains computationally costly at the scales needed for complex plasticity phenomena. Thus, there remains a large gap between what can be studied via first-principles and what is necessary for understanding plasticity problems. Here, we show that newly developed Gaussian Approximation Potential (GAP) for α Fe (bcc, ferromagnetic)²¹ based on machine learning of an extensive DFT data set can contribute to filling this gap.

This GAP potential was been developed in the GAP framework,^22,23 which uses kernel regression and the Smooth Overlap of Atomic Positions (SOAP) to describe the neighbor environment of atoms.²⁴ The GAP potential has been constructed by training it against a large number of atomic environments (~150,000) computed using DFT, including pristine configurations, stacking faults, free surfaces, vacancies, and interstitials. The resulting potential was found to reproduce very accurately DFT vibrational and thermodynamic properties, including equation of state, phonon dispersions, elastic moduli, and thermal expansion.²¹ The potential can be used for system sizes well beyond the capabilities of DFT (up to ~50,000 atoms in this study). Details can be found in ref. ²¹ and are not repeated here. It is also important to remember that while GAP potentials can fail when used to model configurations that are well outside the training set, the underlying framework of GAP potentials includes built-in error indicators to identify such situations. Furthermore, the GAP potential can then be systematically improved by expanding the training set to incorporate additional relevant atomic environments.

Here, we demonstrate that the current generation Fe GAP potential can accurately predict all of the key features associated with the screw dislocations in Fe, relative to DFT. Specifically, the Fe GAP reproduces the DFT-computed core structure, Peierls potential, and slip behavior. Moreover, it shows finite temperature screw dislocation glide by nucleation and propagation of kink-pairs, as envisioned by long-standing theories.^25,26 Importantly, the predicted stress-dependent enthalpy barrier for kink-pair nucleation agrees well, quantitatively and qualitatively, with DFT computations and theoretical models. Overall, this work, together with earlier validations, establishes Fe GAP as the first empirical potential to achieve DFT accuracy for the crucial properties of screw dislocations in Fe, thus enabling future application to many important dislocation/defect problems relevant to the mechanical performance of bcc iron.

The remainder of this paper is organized as follows. Section Results "Dislocation core structure and mobility" contains a benchmark study on the screw dislocation core structure and energetics, showing that the potential developed in ref. ²¹ can be used for simulating kink-pair mechanism. Section Results "Kink-pair nucleation and migration" analyzes the kink-pair nucleation and migration by means of molecular dynamics simulations and nudged-elastic-band computations.²⁷ The main results of this paper are summarized in the Discussion.

Results

Dislocation core structure and mobility

Here, we simulate various aspects of the straight screw dislocation behavior using the Fe GAP via molecular statics (MS) and dynamics (MD) using the LAMMPS package.²⁸ We use a periodic array of dislocations (PAD) configuration (e.g. ref. ²⁹), as described in Methods. The simulation cells have dimensions (l_x = 60a) × (l_y = 2b) × (l_z = 19c), where $a = \sqrt {2/3} a_0$ is the Peierls valleys spacing, $b = \sqrt 3 /2a_0$ the Burgers vector magnitude, $c = a_0/\sqrt 2$ the {101} interplanar spacing, and a₀ the lattice parameter, containing 2,400 atoms. An external stress is applied by assigning forces to the upper and lower Z boundary atoms over a thickness of four atomic layers. For a desired applied stress τ, the forces f and −f on the top and bottom surfaces are

$$f = \frac{{\tau l_xl_y}}{n},$$

(1)

where n = 480 is the number of boundary atoms on the top or bottom layer. Configurations of constant applied stress are computed by relaxing atomic positions while holding the applied boundary forces fixed.

Figure 1 shows the relaxed dislocation core structure at T = 0 K and τ = 0 MPa, projected on to the X–Z plane, along with the associated differential displacement map (DDM). Since screw atomic displacements are along the line of the dislocation, the arrows in the DDM indicate the relative out-of-plane (Y) displacements between the two atoms on either end of the arrow. The size of the arrows is normalized by b/3. Thus, following any closed loop of projected atoms outside the dislocation core yields a total magnitude of displacement of b. In agreement with DFT calculations, the Fe GAP core is compact—a loop around the three central atoms of the core yields magnitude b—and symmetric or non-degenerate, lying in the so-called “easy” core configuration⁶).

Fig. 1

Screw dislocation core as function of an applied stress. Differential displacement maps at the dislocation core are shown for the 〈111〉-zone, as a function of the shear stress τ applied on the free surfaces. Circles represent atomic positions. Colors identify atomic positions in the pristine bcc crystal, for one periodic unit cell (length b) along the dislocation line direction Y, where yellow atoms are taken as reference (Y = 0), red atoms are in position b/3 and gray in position 2b/3. The arrows are normalized by b/3 with a length b/3 scaled to connect neighboring atoms. Arrows represent the displacement of atoms along [1̄11] in the relaxed screw field, computed with respect to the pristine bcc atomic positions

The T = 0 K Peierls stress is computed using molecular statics by incrementally increasing the applied stress τ and relaxing the structure. At the Peierls stress, the dislocation finds no equilibrium configuration and glides steadily through the sample. The computed Peierls stress is τ_P = 2.025 GPa (±5 MPa). Figure 1 also shows the core structure and DDM computed for several applied stresses approaching the Peierls stress. At low stress, there is no visible change of the core structure, consistent with the equilibrium core configuration being a deep energy minimum. As the stress approaches τ_P, the core remains compact. This observation supports DFT-based analyses in bcc Fe and transition metals,^30,31,32 which assume no core “degeneracy”, nor any new intermediate structures as the core shifts under stress, and especially no intermediate split core³² that is a common artifact of most empirical potentials for pure Fe and bcc transition metals.³³

We further compute the Peierls potential, which is the energy change of the straight dislocation line as it moves from the minimum energy configuration towards the next minimum (Peierls valley), at distance a = 0.94b along [121̄], at zero applied stress. Since the core structure distorts as it moves from the local minimum, the Peierls potential is computed using the climbing-image nudged-elastic-band (NEB) method²⁷ (see Methods section), which yields the minimum-energy path (MEP) of the entire system as it moves from an initial state to a final state. The periodic dislocation line length remains 2b, which forces the dislocation to remain straight along the entire MEP.

Figure 2 shows the Fe GAP Peierls potential per 1b dislocation length. The key feature is that the potential has a single hump, and thus no intermediate artificial minimum. This result is consistent with the analysis of the core structures versus applied stress.

Fig. 2

Peierls potential at T = 0 K and zero applied stress. The Peierls potential is shown as computed with the NEB method in the PAD configuration (black line). The Peierls potential computed using DFT in the quadrupolar cell is shown for comparison (red squares), and computations with GAP in the same quadrupolar configuration are also plotted (red line), as reported by ref. ²¹

Also shown in Fig. 2 is the DFT-computed Peierls potential and the corresponding prediction of the Fe GAP model.²¹ Due to size limitations, the DFT computation uses a special periodic quadrupolar unit cell involving two dislocations in which the total dislocation/dislocation interactions, including periodic images, lead to no net glide stress on the dislocations at the equilibrium configuration. The Peierls potential calculation involves the relative motion of the two dislocations and breaks the symmetry, leading to additional energy contributions from the dislocation/dislocation interactions. The result is a larger apparent Peierls potential. The Peierls potential computed using Fe GAP in the same quadrupolar cell is shown in Fig. 2, and the agreement with the DFT is very good, especially considering that the training did not include the Peierls potential but only the gamma surfaces in a 12-atom unit cell. This demonstrates that the true Peierls potential must be computed in a much larger quadrupolar unit cell or by adroit corrections of the direct calculation in small unit cells (see, for example ref. ³⁴).

Other DFT studies report a deduced Peierls potential that is somewhat lower in energy than the results shown in Fig. 2³⁴), and consequently also predict a lower Peierls stress and enthalpy barrier. These differences can be ascribed to different details in the DFT computations (pseudopotential, k-mesh density, convergence tolerances, wave-function cutoff energies) and the type of MEP calculations. Although quantitative differences between various DFT studies are not relevant to the present work, it is useful to point out that the Fe GAP potential was trained on an extensive set of highly converged DFT computations,²¹ making use of a pseudopotential that was shown to reproduce the reference zero-temperature all-electron data.³⁵ Hence, it is expected to faithfully produce behavior consistent with such a reliable computational protocol.

The results above show that the Fe GAP potential developed in ref. ²¹ reproduces the DFT-predicted structure and energetics of straight screw dislocations. This motivates examination of the kink-pair mechanism and computation of the stress-dependent energy barriers for kink-pair nucleation.

Kink-pair nucleation and migration

Having validated the Fe GAP against DFT and conventional understanding for infinite straight dislocations, we now turn to the actual mechanism of glide motion. This problem can only be studied using interatomic potentials, due to the simulation cell sizes required, and so is beyond current full DFT calculations.

Finite-temperature observations of kink-pair nucleation

We first demonstrate that kink-pair nucleation and glide is observed in direct finite-temperature molecular dynamics simulations. Such a study is not influenced by a pre-ordained selection of the final state of the system as imposed in NEB computations described in the next section. We use the same orientation ${\mathrm{X}}||\left[ {12\bar 1} \right]$, ${\mathrm{Y}}||\left[ {\bar 111} \right]$, and ${\mathrm{Z}}||\left[ {101} \right]$ and boundary conditions as specified in Methods. The size of the simulation cell is expanded to contain N = 48,000 atoms with dimensions (l_x = 60a) × (l_y = 40b) × (l_z = 19c) (Fig. 3a). The larger dimensions are necessary to enable nucleation of the correct kink-pair with minimal image stresses due to periodic and free boundaries (Fig. 3b). Stress is applied on the upper and lower Z surfaces as specified in previous Section. MD simulations are performed at finite temperatures. Thus, the lattice constants and cell dimensions are those appropriate to the temperature of interest using the thermal expansion coefficient computed in ref. ²¹

Fig. 3

Dislocation slip at finite temperatures. a Computational cell (black lines), with the atoms out of local bcc symmetry shown (gray), namely the free surfaces (top/bottom) and the dislocation core (in the middle of the simulation cell). b Atoms at the dislocation core during a simulation snapshot, evidencing dislocation glide by kink-pair mechanism. c Atomistic displacement field evidencing slip along the (101) slip plane. Crystallographic visualizations use OVITO⁵⁴ and adaptive Common Neighbor Analysis (CNA),⁵⁵ to label atoms according to local atomic environments

Figure 3 shows the results from one such simulation at T = 200 K. The dislocation begins to glide after 5 ps on the (101) plane of maximum resolved shear stress at τ = 1.3 GPa. The large applied stress is due to the short MD simulation time, and so is not the relevant stress measured in experiments on much longer time scales. The key observation here is that the glide plane is consistent with DFT analyses,^31,32 and the motion is initiated by kink pairs. The kink-pairs nucleate somewhere along the line, and then the kinks glide laterally along the line until, interacting via their periodic images, they attract one another and annihilate, leaving a straight dislocation line that has moved by a along the glide plane. The glide continues by repetition of this unit process. Gliding by kink-pair nucleation is consistent with long-standing theories,^25,26 and kink-pair nucleation and glide have been observed with other potentials.^8,9,10,12,13 This is the first observation of such mechanism with an empirical potential that is quantitatively correct.

Enthalpy barrier versus applied stress

We now compute the stress-dependent enthalpy barrier ΔH^*(τ) for kink-pair nucleation. This barrier controls the rate R of kink-pair nucleation at a given stress τ and temperature T via an Arrhenius law, $R = \nu _0{\mathrm{exp}}\left( {\frac{{{\mathrm{\Delta }}H^ \ast (\tau )}}{{k_BT}}} \right)$, where k_B is Boltzmann’s constant and ν₀ an appropriate attempt frequency. The temperature- and stress-dependent plastic shear strain rate $\dot \gamma$ then follows from Orowan’s law³⁶ as $\dot \gamma \approx \rho b(L/l_{{\rm nuc}})aR$, where ρ is the mobile dislocation density, L/l_nuc is the approximate number of independent kink-pair nucleation sites of length l_nuc available along each dislocation segment of length $L\sim 1/\sqrt \rho$ and aR is the average dislocation velocity. Inverting this relationship leads to the stress as a function of temperature and strain rate.

Due to the work done by the applied field τ acting over the activation area A of the kink-pair nucleus, W = τbA, the enthalpy barrier is a function of the applied stress. To compute the barrier, we again use the climbing image NEB method²⁷ in the 48,000 atom PAD geometry to find the minimum energy path (MEP) followed by the long dislocation line as it moves from an initial state through the transition state to the final state. MEP calculations are done as specified in Methods and Section Results "Dislocation core structure and mobility", with the only difference being that the dislocation length is 40b. Here, fixing the positions of boundary atoms in each replica improves convergence speed, but introduces an error due to the interaction of the non-straight dislocation with the surface. The non-straight dislocation can, however, be viewed as a straight dislocation plus a closed loop of length equal to the kink-pair separation and of width a, and so the image forces decrease fairly rapidly with the simulation-cell size, as ≈(a/(l_z/2))⁻². We have explicitly checked that this boundary effect is negligible by comparing the barriers computed at zero stress with fixed and free boundary atoms in the chosen cell size.

The direct result of the NEB calculation at each applied stress τ is the configurational energy V^C(τ, ξ_i) as a function of the reaction coordinate ξ_i ∈ [0,1] of replica i. This is not the desired enthalpy, as the work done by the applied stress is not included. The NEB enthalpy could be computed atomistically in principle, but this functionality is not available in LAMMPS. To determine the enthalpy from the NEB calculations, the work of the external stress acting on the boundary atoms must be computed for each replica. This work is equal to

$$W^{{\mathrm{ext}}}(\tau ,\xi _i) = {\mathbf{f}}(\tau ) \cdot [{\mathbf{u}}_t(\xi _i) - {\mathbf{u}}_b(\xi _i)],$$

(2)

where u_i (i = b, t) are vectors containing the displacements of all boundary atoms on which forces are applied. As discussed above, this method introduces an error due to the interaction with the boundary, which is minor for the simulation cell size used.

From the replica configurational energy V^C(τ, ξ_i) and the external work W^ext(τ, ξ_i), the replica enthalpy can be calculated^7,8 as

$$H(\tau ,\xi _i) = V^C(\tau ,\xi _i) - W^{ext}(\tau ,\xi _i).$$

(3)

At zero applied stress, this expression coincides with V^C(ξ_i). Figure 4a shows the enthalpy difference ΔH(τ,ξ_i) = H(τ, ξ_i) − H(τ, 0) as a function of ξ_i for a range of applied stress τ from 0 up to the Peierls stress. The path is not of direct importance, but shows the expected behavior. At zero stress there is no strict maximum, because the kink-pair elastic interaction scales logarithmically with kink-pair separation. However, the difference in activation energy between an infinite system and the cell size considered here is negligible. Furthermore, in the finite periodic system, there are kink–kink interactions due to the periodically repeated images of the unit cell and these second interactions cancel due to symmetry when the kink spacing is l_y/2 and hence do not affect the computed enthalpy maximum ΔH^*(0). Therefore, this ΔH^*(0) is also essentially equal to the formation energy for two kinks.

Fig. 4

Energy barrier as function of applied stress. a NEB calculations of the enthalpy difference as a function of the reaction coordinate ξ of the NEB computations for different values of the applied stress. b Barrier height (from the NEB data, markers colored according to a) as a function of applied stress, the line tension (LT) model by Dorn and Rajnak,²⁶ and a Kocks’ law with exponents p = 0.77 and q = 1.3

The main quantity of interest is the enthalpy barrier ΔH^*(τ), which is the peak enthalpy along the minimum enthalpy path. The enthalpy barrier versus applied stress is shown in Fig. 4b. As anticipated, the barrier approaches zero as the stress τ approaches τ_P, and does so in a smooth manner. This general behavior is consistent with previous analyses based on different interatomic potentials for iron.^7,8,9,11,13

The enthalpy barrier for kink-pair nucleation as computed using the Fe GAP can be further compared with the standard models for the kink-pair nucleation mechanism in screw dislocations (Seeger and co-workers,^25,37,38 Dorn and Rajnak,²⁶ Suzuki and co-workers³⁹ and Brunner⁴⁰) recently reviewed in ref. ⁶. We use the early line-tension theory of Dorn and Rajnak, which has some intrinsic simplifications by its neglect of (i) kink–kink interactions, (ii) the dependence of line energy on dislocation character, (iii) the dependence of Peierls potential on the applied stress, and (iv) the asymmetry between left and right kinks.^34,41 However, the model can be calibrated to the computable quantities of the double-kink formation energy, the Peierls stress, and the Peierls barrier (the maximum of the Peierls potential, Fig. 2). Our aim here is to show that our computed ΔH^*(τ) follows the functional trends predicted by this model but typical of all such models.

Dorn and Rajnak²⁶ considered the dislocation as a line defect moving in the Peierls potential V_P(x) under the action of an applied stress. A line energy Γ₀ penalizes changes in length relative to the straight dislocation. Thus, an energy functional for the energy vs dislocation shape can be written down and the saddle point (transition state) configuration can then be determined. The final result of interest here is the enthalpy barrier, given by

$${\mathrm{\Delta }}H^ \ast (\tau ) = 2{\int}_{x_0}^{x_c} \{ [\Gamma _0 + V_P(x)]^2 - [\tau b(x - x_0) + \Gamma _0 + V_P(x_0)]^2\} ^{\frac{1}{2}}{\rm d}x,$$

(4)

where V_P(x) is the Peierls potential per unit dislocation length. Also, x₀ is the equilibrium configuration of the straight screw dislocation at the applied stress τ, such that

$$\tau b = \frac{{\partial V_P(x)}}{{\partial x}}|_{x = x_0},$$

(5)

and x_c is the extremum position of the kink bulge in the critical configuration (transition state) at the applied stress τ, and emerges from equilibrium conditions for the shape of the kink-pair (see ²⁶). The dislocation line energy Γ₀ is related to the kink-pair formation energy ΔH^*(0) via

$${\mathrm{\Delta }}H^ \ast (0) = 2{\mathrm{\Gamma }}_0{\int}_{ - a/2}^{a/2} \sqrt {\left( {\frac{{{\mathrm{\Gamma }}_0 + V_P(x)}}{{{\mathrm{\Gamma }}_0}}} \right)^2 - 1} {\rm d}x.$$

(6)

Here, we fit the Peierls potential V_P(x) as a function of the dislocation core coordinate x along the glide direction (barrier ΔE_P centered at coordinate x = 0) to the standard form²⁶

$$V_P(x) = \frac{{{\mathrm{\Delta }}E_P}}{2}\left[ {1 + \frac{\alpha }{4} + {\mathrm{cos}}\left( {\frac{{2\pi x}}{a}} \right) - \frac{\alpha }{4}{\mathrm{cos}}\left( {\frac{{4\pi x}}{a}} \right)} \right].$$

(7)

We use α = −0.5965 to satisfy the simulated Peierls stress

$$\tau _P = \frac{1}{b}\begin{array}{*{20}{c}} {} \\ {{\mathrm{max}}} \\ {\{ x \in [ - a/2,a/2]\} } \end{array}\left[ {\frac{{\mathrm{d}}}{{{\mathrm{d}}x}}V_P(x)} \right] = 2.025\,{\mathrm{GPa}},$$

(8)

with ΔE_P obtained directly from the simulated Peierls potential (Fig. 2). We then obtain Γ₀ by integration of Eq. (4) with the double-kink energy 2U_k set equal to our computed ΔH^*(0).

The comparison between the computed enthalpy differences ΔH^*(τ) and the classical line tension (LT) model is shown in Fig. 4b. There is close agreement between the classical model and the directly computed activation barrier ΔH^*(τ). The differences reflect the approximate nature of the model, but are minimal. Thus, the results here nicely support the LT model and justify previous applications³⁴ to obtain dislocation mobility laws by fitting to DFT data.

Previous simulations and experiments have also been fit to phenomenological models of the type proposed by Kocks⁴² and widely used for describing thermally activated dislocation behavior in many different circumstances. The general form of the Kocks law is

$${\mathrm{\Delta }}H^ \ast (\tau ) = {\mathrm{\Delta }}H^ \ast (0)\left[ {1 - \left( {\frac{\tau }{{\tau _P}}} \right)^p} \right]^q,$$

(9)

where ΔH^*(0), τ_P, p, and q are parameters that are either fit to experimental/simulation data or are derived from approximate models of the dislocation MEP. Here, we have already computed the values of ΔH^*(0) and τ_P, and so we only need to fit the parameters p and q. The fitted result is shown in Fig. 4b, achieved with p = 0.77 and q = 1.3. These values for (p,q) are close to typical values fitted to experiments on bcc Ta (p = 0.748, q = 1.172⁴³) and bcc W (p = 0.86, q = 1.69⁴⁴). As shown in ref. ⁴⁵ fitting a Kocks’ law to low-temperature experimental data on bcc Fe yields ΔH^*(0) = 0.84 eV, τ_P = 0.363 MPa, p = 0.5 and q = 1. The major difference is not in p and q but in the Peierls stress τ_P, which is ≈5 times smaller than that predicted with either GAP Fe or DFT. This overprediction of the Peierls stress in Fe by DFT and interatomic potentials is well known and is not a specific feature of GAP Fe. Therefore, the GAP Fe potential fits well within the scope of the existing understanding of the activation enthalpy versus stress.

Discussion

We have shown that the Fe GAP potential is suitable for the modeling of screw dislocations in Fe. This is a distinct achievement as the first empirical potential to capture all key features related to screw dislocation glide in bcc Fe, including kink-pair nucleation and propagation at finite stress and temperature. This potential predicts the stress-free compact non-degenerate core structure found in DFT; it predicts the single-hump Peierls potential found in DFT; it does not form common artifact structures such as the split/degenerate core at any applied stresses up to the Peierls stress; and the screw dislocation glides on {110} by the kink-pair mechanism as expected for Fe. Furthermore, the stress-dependence of the enthalpy barrier for kink-pair nucleation is consistent with long-standing theories. These results encourage the use of the Fe GAP potential to study many further issues regarding plasticity phenomena in bcc Fe.

In particular, the deviation between DFT and experimental Peierls stresses at low T has been attributed to nuclear quantum effects⁵⁶ that are absent from DFT and classical simulations. With the present potential, path integral molecular dynamics simulations⁴⁶ could be used to check this hypothesis in the future. Very recent work in this direction has been shown by Freitas et al.⁴⁷. The deviation has also been attributed to the role played by non-screw dislocations when considering expansion of an entire dislocation loop,⁴⁸ and this can also be investigated in detail using the GAP Fe potential.

The Fe GAP potential may also be used to identify the atomic-scale origin of the change of slip mechanism at T ~ 250 K, which was deduced with macroscopic tests in the seminal papers by Brunner and Diehl^49,50 and was more recently observed with in situ TEM.^51,52 Brunner and Diehl proposed that the origin of this phenomenon is the change of elementary slip plane, and hence fitted two models for elementary {110} and {112} slip below and above 250 K, respectively, to reproduce the macroscopic response of experiments. However, despite confirming a change of mechanism, the transition to different elementary slip planes was not observed in the recent in situ TEM by Caillard,⁵¹ who instead interpreted the macroscopic {112} slip traces as being composed of elementary {110} slip events. While GAP Fe has the low-temperature discrepancy in Peierls stress, it can nonetheless probe such slip transitions. GAP Fe can also be employed to predict other parameters that arise in theories, such as the line tension^38,40 and kink-width,^38,39,40 on top of the kink formation enthalpy^38,39,40, which has been computed here (equal to 0.5ΔH^*(0)).

Finally, as noted at the outset, important plasticity phenomena involve dislocation/defect interactions. With the success of GAP Fe in describing what has been the most challenging feature of bcc plasticity—proper behavior of the screw dislocation motion—the study of dislocation interactions with point defects created by irradiation,¹ grain boundaries,² cracks (especially dislocation emission from cracks³), and twinning phenomena⁴ all appear as avenues for future work. We close by reiterating that GAP potentials may fail for new configurations well outside those of the training set but that the GAP formalism enables the identification of such situations and the GAP potential can be systematically improved by expanding the training set to encompass new types of atomic environments arising in such future studies.

Methods

Periodic array of dislocations (PAD) configuration

All atomistic simulations were performed using LAMMPS package.²⁸ To model screw dislocations, we use a periodic array of dislocations (PAD) configuration (e.g. ref. ²⁹). Samples are oriented with ${\mathrm{X||}}\left[ {12\bar 1} \right]$, ${\mathrm{Y}}||\left[ {\bar 111} \right]$, and ${\mathrm{Z||}}\left[ {101} \right]$, with periodic boundary conditions along X and Y, and Z having imposed tractions. We insert a screw dislocation with line direction along Y by imposing a linear displacement u_y = −bx/l_x for 0 < x < l_x on all atoms in the upper half of the simulation cell. Atomic positions are then relaxed by using a combination of the FIRE algorithm⁵³ and relaxation of the cell dimensions until convergence is achieved (force tolerance 10⁻³ eV/atom and stresses σ_XX, σ_XY, and σ_YY <10 MPa). This configuration yields a screw dislocation with periodic boundary conditions along the line direction [1̄11] and glide direction [121̄], and slip plane normal direction [101].

Nudged elastic band (NEB) calculations

In order to perform climbing-image nudged elastic band (NEB) calculations,²⁷ we first compute the relaxed initial state as described in Methods “Periodic array of dislocations (PAD) configuration”. The final state has the same structure as the initial state but shifted by a relative to the initial state. An initial path of intermediate configurations (replicas) is constructed by linearly interpolating the atomic positions between the relaxed initial and final configurations. Replica i is labeled by the reaction coordinate $\xi _i = \left\| {{\mathbf{\xi }}_i} \right\|_2/\left\| {{\mathbf{\xi }}_{{\mathrm{end}}}} \right\|_2$, which is the ratio of the $\ell ^2$-norm $\left\| {{\mathbf{\xi }}_i} \right\|_2 = \sqrt {{\mathbf{\xi }}_i \cdot {\mathbf{\xi }}_i}$ of the 3N-dimensional vector ξ_i = {u_i} containing all atomic displacements of the ith replica (computed with respect to the atomic positions in the initial replica ξ_init) and the vector ξ_end of atomic displacements in the final state. NEB computations are then performed, under the constraint of stress-free boundaries. The force tolerance is set to 10⁻³ eV/Å and the NEB inter-replica spring constant is set to k = 10⁻² eV/Ų. The latter choice does not affect results but optimizes convergence of the calculations.