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

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. A Gaussian Approximation Potential (GAP) can successfully reproduce the structure and motion of screw dislocations in iron. A team led by Francesco Maresca at the EPFL in Lausanne, Switzerland, used the recently developed GAP for iron to simulate aspects of screw dislocation behavior via both molecular statics and molecular dynamics simulations, and validated their results against density-functional theory calculations. The GAP for iron also successfully simulated kink-pair nucleation and screw dislocation glide along the {110} plane, while the stress-dependence of the enthalpy barrier for kink-pair nucleation was consistent with long-standing theories. This potential could be used to identify the atomic-scale origins of many other important plasticity phenomena such as dislocations interacting with radiation damage and cracks in iron and other body centered cubic materials.


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, workhardening, 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. 5 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. 6 Knowing the mechanism, theories exist that can use DFTcomputed inputs. 7 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. 7 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 16 such as the nondegenerate core structure, the single-hump Peierls potential, and slip on {110}. 17 Other features associated with screw dislocations and motion, such as the generalized stacking fault γ surfaces and non-Schmid effects were also investigated. 17 However, the kinkpair 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) 21 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. 24 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. 21 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. 21 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. 21 can be used for simulating kink-pair mechanism. Section Results "Kinkpair nucleation and migration" analyzes the kink-pair nucleation and migration by means of molecular dynamics simulations and nudged-elastic-band computations. 27 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. 28 We use a periodic array of dislocations (PAD) configuration (e.g. ref. 29 ), as described in Methods. The simulation cells have dimensions (l x = 60a) × (l y = 2b) × (l z = 19c), where a ¼ ffiffiffiffiffiffiffi ffi 2=3 p a 0 is the Peierls valleys spacing, b ¼ ffiffi ffi 3 p =2a 0 the Burgers vector magnitude, c ¼ a 0 = ffiffi ffi 2 p the {101} interplanar spacing, and a 0 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 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 6 ).
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-32 which assume no core "degeneracy", nor any new intermediate structures as the core shifts under stress, and especially no intermediate split core 32 that is a common artifact of most empirical potentials for pure Fe and bcc transition metals. 33 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 27 (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.
Also shown in Fig. 2 is the DFT-computed Peierls potential and the corresponding prediction of the Fe GAP model. 21 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. 34 ).
Other DFT studies report a deduced Peierls potential that is somewhat lower in energy than the results shown in Fig. 2 34 ), 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, 21 making use of a pseudopotential that was shown to reproduce the reference zerotemperature all-electron data. 35 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. 21 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 Xjj 121 Â Ã , Yjj 111 Â Ã , and Zjj 101 ½ 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   ). 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 54 and adaptive Common Neighbor Analysis (CNA), 55 to label atoms according to local atomic environments interest using the thermal expansion coefficient computed in ref. 21 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, , where k B is Boltzmann's constant and ν 0 an appropriate attempt frequency. The temperature-and stressdependent plastic shear strain rate _ γ then follows from Orowan's law 36 as _ γ % ρbðL=l nuc ÞaR, where ρ is the mobile dislocation density, L/l nuc is the approximate number of independent kinkpair nucleation sites of length l nuc available along each dislocation segment of length L $ 1= ffiffiffi ρ p 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 27 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)) −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 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 At zero applied stress, this expression coincides with V C (ξ i ). 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.
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, 26 Suzuki and co-workers 39 and Brunner 40 ) recently reviewed in ref. 6 . We use the early linetension 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 26 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 Γ 0 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 (4) where V P (x) is the Peierls potential per unit dislocation length. Also, x 0 is the equilibrium configuration of the straight screw dislocation at the applied stress τ, such that 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 26 ). The dislocation line energy Γ 0 is related to the kink-pair formation energy ΔH * (0) via 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 26 We use α = −0.5965 to satisfy the simulated Peierls stress with ΔE P obtained directly from the simulated Peierls potential (Fig. 2). We then obtain Γ 0 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 34 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 42 and widely used for describing thermally activated dislocation behavior in many different circumstances. The general form of the Kocks law is 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 43 ) and bcc W (p = 0.86, q = 1.69 44 ). As shown in ref. 45 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 56 that are absent from DFT and classical simulations. With the present potential, path integral molecular dynamics simulations 46 could be used to check this hypothesis in the future. Very recent work in this direction has been shown by Freitas et al. 47 . The deviation has also been attributed to the role played by nonscrew dislocations when considering expansion of an entire dislocation loop, 48 and this can also be investigated in detail using the GAP Fe potential.
The Fe GAP potential may also be used to identify the atomicscale 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, 51 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-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, 1 grain boundaries, 2 cracks (especially dislocation emission from cracks 3 ), and twinning phenomena 4 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.

Periodic array of dislocations (PAD) configuration
All atomistic simulations were performed using LAMMPS package. 28 To model screw dislocations, we use a periodic array of dislocations (PAD) configuration (e.g. ref. 29 ). Samples are oriented with Xjj 121 Â Ã , Yjj 111 Â Ã , and Zjj 101 ½ , 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 53 and relaxation of the cell dimensions until convergence is achieved (force tolerance 10 −3 eV/ atom and stresses σ XX , σ XY , and σ YY <10 MPa). This configuration yields a screw dislocation with periodic boundary conditions along the line direction

Nudged elastic band (NEB) calculations
In order to perform climbing-image nudged elastic band (NEB) calculations, 27 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 ξ i ¼ ξ i k k 2 = ξ end k k 2 , which is the ratio of the ' 2 -norm ξ i k k 2 ¼ ffiffiffiffiffiffiffiffiffiffi ffi ξ i Á ξ i p 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 −3 eV/Å and the NEB interreplica spring constant is set to k = 10 −2 eV/Å 2 . The latter choice does not affect results but optimizes convergence of the calculations.

DATA AVAILABILITY
For access to more detailed data than are given in the article please contact the authors.