Temperature and composition dependent screw dislocation mobility in austenitic stainless steels from large-scale molecular dynamics

Extensive molecular dynamics simulations are performed to determine screw dislocation mobility in austenitic Fe0.7NixCr0.3-x stainless steels as a function of temperature ranging from 100 to 1300 K, resolved shear stress from 30 to 140 MPa, and Ni composition from 0.0 to 30.0 at%. These mobility data are fitted to a linear mobility law with a nonzero stress offset, referred to as the threshold stress. We find that both the linear drag coefficient and the threshold stress increase with Ni composition. The drag coefficient increases with temperature, whereas the threshold stress decreases with temperature. Based on these calculations, we determine fitting functions for the linear solute drag coefficient as a function of temperature and composition. The mobility laws determined in this study may serve to inform dislocation dynamics simulations pertinent to dislocation network evolution at elevated temperatures for a wide composition range of austenitic stainless steels.


INTRODUCTION
The use of austenitic stainless steels such as 304 and 316 is ubiquitous across applications that demand strength and corrosion resistance at moderate cost. Advanced austenitic steels are undeniably among the most important materials of the modern age, spanning automotive, aerospace, nuclear, and consumer household applications. Recently, advances in nonequilibrium additive manufacturing methods such as selective laser melting (SLM) have revitalized scientific inquiry into processing and microstructural properties of such steels. Specifically, the formation of complex 3D dislocation networks and chemical ordering has been observed, leading to novel deformation behavior without the traditional tradeoff between strength and ductility [1][2][3][4] . The phenomena that affect material responses and properties span multiple length scales. The modeling of such defect structures is an inherently multiscale problem 5,6 , and balancing accuracy and applicability becomes increasingly important. The behavior of individual dislocations, microstructure-scale mechanisms, and collective dislocation interactions must be bridged to enable predictive simulation of dislocation plasticity. Discrete dislocation dynamics (DDD) is one such method that seeks to describe mesoscale phenomena by evolving microstructures using a set of physics-based heuristic rules [6][7][8] . Continued advancements in the implementation of such rules has equipped DDD with the ability to describe mechanisms such as dislocation glide, climb, and cross slip, which contribute to development of threedimensional dislocation structures 9,10 . A fundamental driver in all such extensions is the mapping of the local resolved shear stress to the dislocation glide velocity on specific slip systems, described by the mobility law [11][12][13][14][15] . While difficult to ascertain experimentally, mobilities have been readily computed via molecular dynamics (MD) simulations for pure metals and binary alloys for both face-centered-cubic (fcc) and body-centered-cubic (bcc) crystals [12][13][14][15][16][17][18] . To date, the lack of mobility data for more complex alloys has limited the application of DDD simulations to the SLM austenitic stainless steels of interest, with some authors 19,20 opting to take single bulk handbook values for the mobility despite clear evidence of temperature and composition dependence 12,16,21 . Furthermore, the temperature ranges simulated in typical mobility calculations must be extended to capture the high temperatures induced by the laser path traversal and subsequent cooling during SLM. Moreover, austenitic stainless steels are often employed in elevated temperature applications and these ranges are simultaneously captured. Lastly, an examination of composition effects on dislocation mobility is key to understanding the effects of alloying on the strengthening behavior. In this work, we conduct an extensive MD study on screw dislocation mobility in the austenitic Fe 0.7 Ni x Cr 1−x stainless steels in an extended temperature and composition regime.
The dislocation mobility can be generally divided into thermally activated 17,21 regime, a linear drag-dominated 16,22 regime which at higher strain rates represents phonon drag, and an asymptotic shear wave-limited relativistic regime 12,23,24 . In pure fcc metals, where the Peierls stress (i.e., minimum stress to move the dislocation at 0 K) is negligible, only the phonon and relativistic regimes are observed. Mobility functions relate the resolved shear stress τ to the dislocation velocity v, and typically reflect a piecewise phenomenological form, i.e., Here, fitting parameters A and D are linear and nonlinear drag coefficients, respectively, b is the Burgers vector length, and v 0 a minimum critical velocity for radiative dispersion. It is also common to define a friction coefficient that incorporates temperature effects as AT = B(T), though this term is often interchangeable with drag coefficient 17,25,26 . The power law exponent of 3/2, as originally suggested by Eshelby for the dispersive term 23 , is commonly cited in the literature; however, recent studies have sought to render more accurate fits as a function of character angle and stress state 13,14,27 . Moreover, various interpolations between screw and edge character drag coefficients are typically used in the absence of this data during conventional DDD simulations 28 . In the low stress regime at zero temperature, a threshold stress may be required to initiate dislocation glide in the fully linear drag-dominated regime. For a pure metal this is governed by the Peierls stress. In random solid solutions, the presence of solute-dislocation and solute-solute interactions 21 induces elastic misfit strain fields 29,30 , presenting non-negligible pinning barriers 16 that must be overcome to initiate dislocation glide. This threshold stress controls the material's resistance to dislocation glide at low temperatures, where the kinetics of thermally activated obstacle bypass is slow. Recently, efforts to predict this in the context of concentrated solid solutions have been proposed 30 and will be discussed later in this work. However, most DDD simulations in the literature consider only the regions described by Eq. (1), applying the linear drag fit through the regime of lower effective stress range 14,19,28,31,32 . On the other hand, the threshold for the onset of dislocation glide becomes increasingly negligible at elevated temperatures 33 at which thermally-assisted dislocation migration is more pronounced. In this work, we focus primarily on the linear drag region while characterizing variations in the threshold stress. Furthermore, for the range of applied resolved shear stress explored here we omit dispersive (high velocity) effects 34 .
In order to describe temperature and composition effects on linear drag induced by the solute fluctuation in a binary alloy, Marian and Caro 12 proposed an additional term in the drag coefficient A, i.e., where x is the solute concentration. This suggests a linear dependence on composition. However, this form was derived and validated against a binary alloy within a limited temperature range; the presence of multiple solute atoms may require additional terms or an alternate form altogether, particularly if their relative misfit volumes differ 30 . In most previous work seeking to define DDD mobility functions, explicit values for drag coefficient have been calculated and used at a few select temperatures and compositions of interest. To the best of our knowledge, the current work offers the most comprehensive analysis to date across both temperature and composition variables in the Fe 0.7 Ni x Cr 0.3−x austenitic stainless-steel alloy system.  Fig. 1a for temperature = 100 K and Fig. 1b for temperature = 500 K. Two linear fits to each set of MD data points are presented: the green line indicates a linear fit excluding velocities below 0.1 Å/ps, while the red line indicates the linear fit through all data points. For subsequent analysis we fit to the truncated mobility data defined here. The R-squared value (inset, Fig. 1) of these fits decreases with increasing temperature. This can be attributed to the nonnegligible threshold stress from the thermally-activated regime to the onset of drag-dominated dislocation glide, as discussed in "Methods" section. Significant threshold stresses at low temperature are observed across all compositions. When temperature is increased to 500 K, it is observed that this threshold stress is diminished, and a linear approximation through the entire stress range is increasingly valid at higher temperatures. We do not enforce a fit through the origin; the regression is allowed to compute a nonzero stress intercept that would correspond to the transition to the thermally-activated regime at lower stress levels. A detailed discussion regarding the use of least squares regression in determining drag coefficient can be found in the ref. 22 .

RESULTS AND DISCUSSION
Dislocation mobility curves Average dislocation velocities (over the five replicas) are calculated as a function of shear stress for all compositions and temperatures simulated, and results are summarized in Fig. 2a-g. The more significant nonlinearity at low temperatures (100 and 300 K) and the nonzero stresses for the onset of dislocation motion in the linear solute drag regime seen in Fig. 1a, b, can be observed across all compositions. It is also apparent that the slopes of the mobility curves decrease with increasing temperature at all compositions, consistent with drag-controlled dynamic behavior. In other words, the drag coefficient increases with temperature once the initial threshold stress is overcome; this is discussed further in the following sections.
Qualitatively, Fig. 2 shows no evidence of a plateaued relativistic or dispersive regime for the stress range tested. Indeed, the critical velocity v 0 for crossover into this regime is comparable to the lowest shear wave speed in the solid 12 , and was determined for various other systems to be at least double that of our highest observed velocity of~8 Å/ps 13,16,22,32,35 . Taken together, these two observations indicate that a linear mobility law with nonzero stress at zero velocity suffices in the viscous damping region.
To examine solute effects, we plot the average velocity as a function of Ni composition for different stresses and temperatures, shown in Fig. 3. It is evident that the average velocity decreases monotonically with increasing Ni content, suggesting a stronger pinning contribution from Ni relative to Cr. Note, we also considered binary Fe 1−x Ni x alloys and a similar decrease in velocity was observed as the alloy concentration increased (see Supplementary Fig. 3).
Local environment effect on dislocation mobility The dislocation mobility in random alloys is controlled primarily by solute-dislocation interactions. Local environment heterogeneity resulting from fluctuations in swelling volume and elastic constant differences are important to capture in executing MD simulations of Composition dependent dislocation velocity. Averaged velocity in the ternary alloy as a function of Ni composition. Lines plotted here serve as a visual guide to demarcate data and do not reflect any fitting procedure.
f) e) Fig. 2 Dislocation mobility curves. Dislocation velocity as a function of resolved shear stress for different Ni compositions at temperatures a 100 K, b 300 K, c 500 K, d 700 K, e 900 K, f 1100 K, and g 1300 K. Lines plotted here serve as a visual guide to demarcate data and do not reflect any fitting procedure.
this nature. To expand on this, average atom simulations are executed for a fixed composition and temperature of x Ni = 0.15, T = 300 K for comparison. The simulation cell and parameters follow those outlined in "Methods" section. The dislocation position is plotted as a function of timestep for a single simulated stress τ = 60 MPa in Fig. 4a, and the full mobility curve for these parameters is shown in Fig. 4b. It is evident here that the stress required to initiate dislocation motion is much lower for the average atom results, resembling a pure FCC crystal. In addition, the glide velocity of the dislocation is also significantly higher for a given applied shear stress when using a purely average atom representation. In fact, many solute-strengthening theories like the one explored later in this work, among others 21,24,29 , derive from the interaction between solute and dislocation stress fields during glide. While an average-atom potential may very well capture the static energetics and bulk properties of alloys, these direct comparisons show it does not adequately capture glide in random alloys. This fact was pointed out in 36 , though a mixed average atom potential wherein true solutes interact with average atom types may address this. A detailed comparison of the lattice distortion effect stemming from random solute distributions for a different system can be found in the ref. 37 .
where c is introduced to represent the intercept computed by the fit. Here, by applying the regression through the truncated mobility data as described in Fig. 1, and the drag coefficient at each simulated temperature and composition can be extracted. Figure 5a, b show the computed drag coefficients as a function of temperature and Ni content and vice versa. The drag coefficient increases significantly with increasing temperature and Ni composition. Figure 5a shows that the dependence of the drag coefficient on temperature can be reasonably well approximated as linear. The fitting parameters for Fig. 5a are reported in Table 1.  Table 1. Table 1. Parameters for drag coefficient as a function of temperature:

Threshold stress
In solid solutions, dislocations are pinned at obstacles and must be overcome by the applied stress in order to initiate glide. This threshold stress τ th is sometimes referred to as the static friction stress 12,16,30 ; we define this threshold as the resolved shear stress at which the dislocation velocity exceeds 0.25 Å/ps, which is approximately the smallest velocity which could be obtained from the MD data 38 . As emphasized previously, the mobility curve does not remain linear for the entire stress and temperature range, in particular at the lower temperatures considered; therefore, the stress corresponding to 0.25 Å/ps is determined by interpolation of the MD data. At elevated temperature, where the mobility curve is linear, we estimate the threshold stress as the x-axis intercept value from the linear fit to the MD data. A comparison of intercept and interpolated threshold stress is shown as a function of temperature in Fig. 6a-g for different Ni compositions. The friction stresses are fit to an exponential function of temperature of the form τ th = A 0 exp(−δT); the fits are shown in red in Fig. 7. Note that the interpolated stresses are used for the fit where possible. We can thereby estimate the friction stresses at 0 K, τ th0 , by extrapolating the exponential fit. Fitting parameters for each composition are listed in Table 2.
Solute effect on drag coefficient In simulations of complex alloy systems, dislocation glide is sensitive to the random seed used to generate the alloys 21 . To address this, the dislocation velocities were averaged and large simulation cells were used. Additionally, our calculations (Supplementary Fig. 2) and previous work 38 indicate that there exists a critical dislocation line length below, which dislocation mobility is sensitive to the simulation cell size due to the pinning-depinning characteristics of dislocation motion in alloys, especially at lower stresses (see refs. 21,35,38 for a more detailed discussion of this phenomena). In this work, the dislocation line length is~400 nm, almost ten times longer than that used in ref. 21 , and five times larger than that used in ref. 12 . Together, these improvements lend confidence to the composition-dependent results acquired here. Moreover, the drag coefficient dependence on solute concentration has been clearly identified by various MD simulations in Interpolated stress Intercept stress Exponential fit   binary alloys 12,16,21 . In previous work on pure elements, it has been observed that the drag coefficient is nearly independent of temperature, when calculated as A = B/T in Eq. (1) 16 ; however, when our drag coefficient calculations are normalized by temperature this is not observed, particularly at lower temperatures. A similar behavior was observed in 38 , where a power-law dependence of the drag coefficient on temperature was necessary to fit the MD data. A possible explanation for this is the additional variation in energetic barriers introduced by solute atom interactions with the dislocation. This temperature dependence can become more complex with an increasing number of components, since multiple solute types can have interacting stress fields. The drag coefficient expression, Eq. (2), derived by Marian and Caro 13 was initially based on narrow ranges of concentration and temperature for binary systems. In this equation, the composition and temperature effects are coupled into a single multiplicative term, so that the intercept parameter A * reduces to a temperature-independent and compositionindependent value for the pure system. In contrast, we cannot fit our data to Eq. (2) because the system explored here is a ternary alloy, resulting in multiple values for the intercept; therefore, we point to Table 1 for a more easily interpreted form of the mobility function. We may still estimate the derivative dB/dx of the drag coefficient with respect to concentration x at each temperature and obtain β(T), as defined in Eq. (2), from the linear regressions shown in Fig. 5b. Note that since we are simultaneously varying x Ni and x Cr , dB/dx = dB/dx Ni = −dB/dx Cr and our fit to β(T) is only strictly valid for the case of x Ni + x Cr = 0.3 considered here. The result is shown with its linear fit in Fig. 7.
The above analysis suggests that Eq.
(2) is difficult to apply to ternary alloy systems and other complex alloys (e.g., high-entropy alloys). Figure 5 reveals that the drag coefficient remains fairly linear as a function of both variables; therefore, we explore a multiple regression approach: Here, ω and γ are linear fitting constants. Equation (3) decouples the temperature and composition dependence but enables a single compact drag coefficient function to be implemented into DDD simulations for alloys with a wide range of temperature and composition. Its linear nature does not conflict with the existing theory, regarding temperature and composition, on dislocation mobility.
Parameters of Eq. (3) are fit to our MD data and are listed in Table 3. Using the fitted parameters, the estimated drag coefficient is shown as a function of temperature and composition in Fig. 8a. The absolute percent error of the fitted expression is shown in Fig. 8b. We can see that the functional form in Eq. (3) results in a good fit for the entire domain. (maximum error of 10% at low temperature and an average error below 3%). This confirms that the temperature and composition effects are not strongly nonlinearly coupled in the ranges investigated and separation of these two variables is a good approximation.
Theoretical predictions of zero temperature threshold stress The 0 K friction stress calculated above is the purely static contribution of the solute atoms to material strengthening in the absence of kinetic contributions at finite temperature. Conceptually, it is analogous to the zero-temperature flow stress τ th0 recently derived by Leyson et al. 39 for solid solutions and applied to high entropy alloys by Varvenne et al. 30 . Based on the purely elastic interactions between the solute and dislocation's hydrostatic stress field, Varvenne et al. derived an analytical expression for τ th0 as : (4) Here, x n is the concentration of solute species n, α is the dislocation line tension parameter, μ is the alloy shear modulus, ν is the Poisson's ratio, f 1 (w c ) is a factor describing the bowing geometry, b is the magnitude of the composition-weighted Burgers vector component in the glide (edge) direction, and ΔV n and σ ΔVn are the mean misfit swelling volume and associated standard deviation, respectively. f 1 is referred to as the minimized core coefficient and was derived by Varvenne et al. 30 to be a constant 0.35 for the fcc edge dislocation case in Fe 0.70 Ni x Cr 0.3-x . In general, the partial spacing is reduced for pure screw dislocations 30 , and therefore the value of f 1 must be recalculated for screw dislocations. Following the approach of Varvenne et al. 30 , we recalculate f 1 = 0.059 (see Supplementary Figs. 4, 5). The line tension parameter is α = 0.492 for screw dislocations based on MD calculations of screw and edge line tension in Al 40 . Elastic constants μ and ν are determined by molecular statics at 0 K for each solute concentration in the ternary alloy using five random alloy configurations for each concentration with 4000 atoms in the simulation cell. We calculate the "isotropic average" elastic constants by rotating the Voigt elastic tensor through all orientations normal to the dislocation Burgers vector ( 110 Â Ã direction) 41 . To illustrate, Fig. 9a shows Poisson's ratio as a function of the rotation angle and Fig. 9b shows the averaged Poisson's ratio and shear modulus μ as a function of Ni composition. Table 3. Parameters for drag coefficient as a function of temperature and composition: B(T, x) = B 0 + ωT + γx. In order to calculate τ th0 using Eq. (4), ΔV n and σ ΔVn are needed. Following the previous approach 38 , molecular statics simulations with 868 atoms are performed to calculate ΔV n and σ ΔVn associated with substitution of a single solute atom for every lattice site in the system. We loop over all lattice sites in a random ternary solid solution at each target composition for both Cr and Ni solutes. When the site is occupied by Fe, it is switched to the desired solute atom. When the site is occupied by the solute atom, it is switched to Fe. The average volume change ΔV n due to the insertion of the solute atom to a random Fe 0.7 Ni x Cr 0.3−x alloy can thus be determined. The resulting ΔV Cr , ΔV Ni , ΔV Fe values and their standard deviations are shown in Fig. 10 as a function of x Ni . Using ΔV n and σ ΔVn , the 0 K friction stress is calculated using the theoretical expression Eq. (4).
The results of τ th0 obtained from the theory and MD are compared in Fig. 11. The molecular statics results indicate large standard deviation σ ΔVn , in some cases exceeding the computed mean value ΔV n . This observation is unsurprising due to the highly variable local solute environment in the alloy system. The friction stress predictions shown in ref. 30 neglected the influence of σ ΔVn for each solute, whereas our computations here explicitly include it. We similarly calculate μ and ν at each composition rather than taking an averaged value, and the elastic constant terms introduce another degree of composition dependence that is reflected in the analytical model results (Fig. 11).
We note that the predictions presented in Varvenne et al. 30 below a temperature of 200 K also deviate from the available experimental results by over 50 MPa for the alloys investigated, and Varvenne et al. 30 noted that the predictions become increasingly dependent on the line tension parameter α at low temperatures. In addition, while we rederived the core parameters for the dissociated screw dislocation, the edge components cancel beyond the core region, and the long-range stress fields for the screw dislocation do not interact with the purely hydrostatic solute misfit stress. The Varvenne et al. theory accounts for solute interactions at both short and long ranges; however, so we do not expect this to affect the applicability of the theory for screw dislocations.
The final point is that this theory was developed and validated against equiatomic multicomponent alloys, whereas the present work considers a relatively dilute solid solution in comparison. As such, the use of averaged values in the analytical expression may be sufficient for high entropy alloys; however, in our case we notice that the Varvenne et al. 30 expression predicts a stronger dependence on composition which is not reflected in the MD results. One interesting feature that is captured by both models is the non-monotonic behavior at x Ni = 0.25, which emphasizes a) b)  again that the fluctuations in solute interaction with the dilatational stress field of the screw dislocation is a major contributor to strengthening. Extensive MD simulations carried out in this work have allowed us to obtain reliable mobility data for the Fe 0.7 Ni x Cr 0.3−x alloy. We find the following: 1. The dislocation drag coefficient follows a consistent trend with temperature and composition in the linear solute drag regime. In particular, Ni has a more pronounced dragging effect than Cr on screw dislocations. This effect is slightly stronger at higher temperatures (Fig. 4). 2. The conventional drag coefficient B(T) has been computed as a function of temperature at each composition for ternary Fe 0.7 Ni x Cr 0.3−x alloys representative of austenitic stainless steels. This is useful to inform DDD simulations at a fixed alloy composition. 3. We have proposed a bivariate linear fit of the drag coefficient applicable to temperature, stress, and composition range. The error of the fit is low for applied stress between 10 and 140 MPa, suggesting weak nonlinear coupling between composition and temperature effects in this range. 4. Our results are compared to a recently developed analytical model for strengthening in concentrated solid solutions which we extend to the screw dislocation case, suggesting predictive limitations for low temperature behavior, notably the 0 K friction stress.
As a direct outcome of this work, the functional form for the drag coefficient as a function of temperature and solute composition in austenitic stainless steels may be implemented directly in mesoscale models. This is an important advancement in the development of a physically informed mobility model for enabling accurate DDD simulations of 3XX austenitic stainless-steel. In particular, it is applicable to elevated temperatures and a range of composition typical in applications of interest or associated with composition gradients resulting from additive manufacturing processes. The high cost of simulating such large systems and parameter spaces at full atomic resolution may motivate the use of coarse-graining techniques such as the concurrent-atomistic-continuum (CAC) method 42 or other multiscale methods 5 in future studies, particularly for considering extended defect-defect or defectinterface reactions.

Interatomic potential
Herein we employ a recently developed Fe-Ni-Cr embedded atom method (EAM) potential 43 as a surrogate for austenitic stainless steel. This potential was selected because: (1) it enables stable MD simulations of the fcc structure; (2) it captures the correct unstable stacking and stable stacking fault energies (SFE); (3) it accurately reproduces elastic constants; (4) it gives reasonable trends of various energies and volumes for a range of compositions; and (5) it predicts that solutes remain in solution for the composition range of interest here. The latter two points are particularly relevant when considering the overall effect of alloying. Together, these five features are pertinent to studying the mechanical behavior of stainless steels. In the context of dislocation behavior, the correct SFE trend with composition is especially important, as the SFE controls equilibrium Shockley partial spacing. The potential selected matches closely the trend observed in both DFT and experiments [44][45][46] (see Supplementary Figs. 6, 7). Phonon characteristics for this potential can be found in the ref. 38 and are important in controlling phonon drag behavior. One thing to note is that that while varying the composition in Fe 0.7 Ni x Cr 0.3−x , the alloy is not strictly austenitic at x Ni = 0.0 but remains fcc at all temperatures considered here. This does not affect the interpretation of composition-dependent results obtained. Average-atom EAM potentials 36 have also recently been proposed as an efficient way to represent random alloy characteristics.
We generate average-atom potentials for each target composition studied here to compare with the true random alloy results reported.

Molecular dynamics simulation
Our simulation cells embedded with a single screw dislocation contain 11,108,370 atoms with dimensions of approximately 4000 Å in x, 160 Å in y, and 200 Å in z. This large dimension in the dislocation line direction allows for additional sampling of a wide range of solute environments, giving a more accurate dislocation line response and size-independent velocities 38 . Using different random number seeds, five replica simulations were performed for each stress and temperature condition to reduce any effect of the particular initialization method.
As shown in Fig. 12, a displacement field is applied to the atoms to create a perfect screw dislocation with both its line direction and Burgers vector b * ¼ ½110a=2 parallel to x. The initially perfect dislocation is allowed to relax, resulting in its dissociation into two Shockley partials on the expected (111) slip plane.
Our simulations employ periodic boundary conditions in x and z, and a free boundary condition in y. Note that for this fcc geometry, the stacking period in z is 6(112) planes as in the ABCDEFABCDEF… sequence. Our nominal size (dislocation line direction) is set as large as possible to circumvent the dislocation length-dependent regime for most simulated conditions (stress and temperature) 38 . Weinberger et al. 22 have previously investigated the effect of nonperiodic line length; however, this introduces the effects of pinning at the free surfaces, which we avoid here by using fully periodic boundaries in both the glide(z) and line(x) directions. The effects of y and z (nonline direction) dimensions were not explicitly studied previously 38 and are included in Supplementary Fig. 1; we found that the y and z dimensions do not affect the velocity in an appreciable way. Moreover, no cross slip along the screw dislocation was observed for any of the configurations simulated.
During the MD simulations, the upper and lower 2(111) planes of surface atoms are subject to constant opposing forces parallel to the Burgers vector (i.e., parallel to x) to obtain the intended resolved shear stress is τ. The initial velocities of atoms in these two surface layers are set to zero, and initial velocities of other atoms are assigned based on a Boltzmann distribution and the simulated temperature. Atom trajectories are then solved using an NVE (constant volume and energy) ensemble for the surface atoms and a 1 atm pressure NPT (constant pressure and temperature) ensemble for the other atoms. All MD simulations are performed for a total of 4 ns using LAMMPS 47 , and dislocation position as a function of time is calculated from total relative displacement of the crystal above and below the dislocation glide plane following the previous approach 38 . We employ standard thermostat/barostat damping parameters of 100 * dt and 1000 * dt, respectively, with a timestep size of dt = 0.004 ps. To ensure that the results were not sensitive to the choice of damping parameter, mobility calculations were also run for a fixed configuration with three different values for each. There was negligible variation between the velocities calculated. Fig. 12 MD simulation cell geometry. MD simulation cell geometry for dislocation mobility simulations. The initial position of the inserted screw dislocation is indicated by the black arrow. Atoms in green are fcc coordination while white atoms indicate the free surface in the y direction.