PHYSICAL REVIEW B 87, 144106 (2013)
Prediction of the kink-pair formation enthalpy on screw dislocations in -iron by a line tension model parametrized on empirical potentials and first-principles calculations
Laurent Proville,1 Lisa Ventelon,1 and David Rodney2
1
CEA, DEN, Service de Recherches de M´ etallurgie Physique, F-91191 Gif-sur-Yvette, France 2 SIMAP-GPM2, Grenoble INP, CNRS/UJF, F-38402 Saint-Martin-d'H` eres, France (Received 14 January 2013; published 23 April 2013)
The thermally activated glide of screw dislocations in bcc crystals proceeds through the formation of kinks by pairs, a mechanism named after Peierls for his early work on dislocation theory. A method is proposed to compute the dislocation kink-pair formation enthalpy from density functional theory (DFT) calculations. This method consists of properly adjusting the parameters of a one-dimensional line tension model from atomistic calculations performed in small simulation cells. This model is applied to bcc iron to determine the kink-pair formation enthalpy at different applied stresses from DFT calculations. DOI: 10.1103/PhysRevB.87.144106 PACS number(s): 62.20.F-
I. INTRODUCTION
The temperature-dependent plasticity of bcc crystals is controlled by the thermally activated glide of a0 /2 111 {110} screw dislocations.1,2 These dislocations are the purpose of great computational efforts (see Refs. 316 for examples), the main issue being deriving the dislocation mobility law from the atomic scale.17 Among the different computational studies, calculations based on the density functional theory (DFT) have established some important features that were previously matters of debate in transition metals, such as the nondegenerate dislocation core structure,6 the {110} glide plane, and the single-hump Peierls barrier.10 However, the main drawback of DFT-based methods applied to extended defects such as dislocations is the limited number of atoms that can be considered in the simulation cell, i.e., typically a few hundred. Such numbers are sufficient to study straight dislocations perfectly aligned with their b = a0 /2 111 Burgers vector because the translational invariance along the dislocation line allows us to reduce with periodic boundary conditions the cell size in this direction to a slab of length b = |b| = a0 3/2.912 At low but finite temperatures, however,18 the motion of screw dislocations proceeds through the formation of kink pairs, thereby breaking the translational invariance along the dislocation. Simulating this process at the atomic scale thus requires us to consider a long three-dimensional dislocation with kinks, which is out of the reach of current DFT calculations. The increase in dislocation length induced when the dislocation locally bows out to form a kink pair is associated with a restoring force, called line tension (LT),1922 in reference to the standard picture where the dislocation is viewed as an elastic string. This analogy is not exact in particular because it is based on a local approximation which is not adapted to dislocations in the general case. The reason is that dislocations interact at long range with themselves and other dislocations. However, we show here that the line tension approximation is applicable to kink pairs. On the other hand, our work underlines the fact that the classical expression derived from linear elasticity23,24 cannot be used for kink pairs. Indeed, this expression depends on the dislocation configuration considered with a geometric term ln( /r0 ), where is the length of the bowing segment and r0 is the inner cutoff distance, typical of linear elasticity,
1098-0121/2013/87(14)/144106(8)
of the order of the Burgers vector. When a kink pair nucleates, the distance between kinks, which sets the size of the bowing segment, is initially on the order of the cutoff distance, so that the uncertainties on the line tension are in the same range as the physical effect to be captured. Thus, as for other dislocation features mentioned above, the dislocation line tension requires us to account for the details of the atomic interactions. In this study, we present a method to compute the dislocation line tension in atomic-scale simulations in the case of screw dislocations in bcc crystals. The line tension calculations are performed in a simulation cell small enough to be accessible to DFT calculations, performed using the SIESTA code.25 The line tension is then used as a parameter in a line tension model19,26 to predict using DFT calculations, the kink-pair formation enthalpy on screw dislocations as a function of applied stress, a calculation which is beyond the reach of direct DFT calculations. In order to check the validity of this method, the same calculations were performed using three different semiempirical interatomic potentials for iron that are all well adapted to simulate screw dislocations. Although less accurate than DFT,27 the advantage of these potentials is to allow for a direct calculation of the kink-pair formation enthalpy on three-dimensional screw dislocations that we compared to the LT predictions.
II. ATOMISTIC MODELS AND SIMULATION CELLS A. Embedded-atom-method potentials
The three embedded-atom-method (EAM) potentials employed in this work to model bcc iron were presented in previous studies. The first potential, denoted here as EAM 1, is the well-known and widely used potential developed by Mendelev et al. in Ref. 28. In the latter reference, this potential was denoted number 2. The second potential, denoted EAM 2, is the potential named MPG20 by Gordon et al. in Ref. 29. This potential was developed to remedy the main shortcoming of EAM 1, which is a pronounced metastable split state for the screw dislocation in between Peierls valleys. The third potential, denoted EAM 3, was recently proposed in Ref. 16 to correspond more closely to DFT calculations, in particular regarding the shape and height of the Peierls barrier. All three potentials predict a nondegenerate core configuration for the
©2013 American Physical Society
144106-1
LAURENT PROVILLE, LISA VENTELON, AND DAVID RODNEY
PHYSICAL REVIEW B 87, 144106 (2013)
screw dislocation, a {110} glide plane at low temperature and a Peierls stress of about 1 GPa, in agreement with DFT calculations.6,10,14,30
B. DFT calculations
The present first-principles electronic structure calculations were performed within the DFT framework using localized basis sets as implemented in the SIESTA code.25 The SIESTA pseudopotential and basis set for iron were previously validated by stacking fault energy and dislocation calculations.10,11 The 2 × 2 × 16 and 2 × 2 × 8 shifted k -point grids were used for the 1b and 2b cells, respectively. Residual forces after ° The Hermite-Gauss relaxation were smaller than 0.01 eV/A. scheme to broaden the electronic density of states was used with a smearing of 0.3 eV. The calculations in iron are spin polarized (ferromagnetic iron), and the generalized gradient approximation (GGA) was used.
C. Simulation cells
quadrupolar geometry does not allow us to apply constant stresses in nudged-elastic-band (NEB) calculations. As a result, in the following, all DFT calculations were performed with no applied stress. Keeping in mind that DFT simulations were performed with xz = 0, we use the same notations as those introduced previously for the EAM case.
III. LINE-TENSION MODEL A. Formulation and approximations
The simulation cell is orientated such that the screw ¯ plane, the dislocation glide plane is a horizontal Z = [110] dislocation line is along the X = [111] direction, and the ¯ 12] ¯ axis. Depending on the glide direction is along the Y = [1 energetic model, EAM or DFT, different cell geometries and boundary conditions were used in order to minimize finite-size effects in particular in DFT calculations where the number of atoms is very limited. In the case of EAM simulations, we considered a single dislocation with periodic boundary conditions in the X and Y directions and free boundary conditions in the Z direction.31 ¯ The upper and lower free surfaces are [101] planes of the crystal, denoted S + and S - , respectively. To reproduce the effect of an applied shear stress xz ,31,32 external forces are added to each atom in S + and S - , fx+ = xz ab and fx- = -xz ab , respectively (ab is the surface per atom in S ± , with a = a0 2/3 being the distance between Peierls valleys). The total number of atoms is denoted Nat . The total energy Vat ({ri }) is the sum of the interatomic potential energies, with {ri } being the 3Nat vector of atomic positions. The enthalpy of an atomic configuration is expressed as Hat ({ri },xz ) = Vat ({ri }) - xz ab
j S +
Modeling a dislocation through a one-dimensional (1D) elastic line on substrate is a standard approach in dislocation theory, presented as the 1D Frenkel-Kontorova model3335 or as the string model36 or, more recently, as the line-on-substrate model.37 It is also often simply named the line-tension model.19,26 The screw dislocation is then represented by a single function y (x ), which is defined as the position y in the dislocation glide plane of the dislocation segment situated at the x coordinate along the dislocation line. If we assume that the derivatives of the function y (x ) are small and furthermore that the Peierls barrier is small compared to the dislocation self-energy, we can write a linearized version of the LT model,3739 where the dislocation enthalpy is expressed as HLT (y,xz ) = dx HP (y (x ),xz ) + T 2 dy dx
2
.
(2)
xj -
j S -
xj , (1)
where xj S + (xj S - ) is the X coordinate of atom j in S + (S - ). In the following, enthalpies at a given applied stress are computed with respect to a reference configuration, {r 0 }, where the dislocation is straight and relaxed under the action of both the interatomic potential and the applied stress, i.e., the minimum of Eq. (1). In the case of DFT calculations, we considered a quadrupolar geometry, where a dislocation dipole was introduced in a simulation cell with periodic boundary conditions in all three directions.10 The reason for this different geometry is that a quadrupolar arrangement with periodic boundary conditions allows us to correct for finite-size effects and can thus be used with very small numbers of atoms, smaller than with the single-dislocation geometry employed in the EAM case. The drawback is that, in the absence of free surfaces, the
Here, T is the line tension, which in this simple model is assumed constant, and HP is the substrate enthalpy. If the dislocation Peierls potential, denoted VP , does not depend on the applied stress, the substrate enthalpy can be written as a linear function of the applied stress: HP (y (x ),xz ) = VP (y (x )) - xz by (x ). However, in some cases, atomistic calculations showed that the Peierls potential depends on the applied stress,38,40 reflecting either an evolution of the core structure under stress or a stress dependence of the dislocation pathway between Peierls valleys. Therefore, we prefer to keep the Peierls potential and the work of applied stress together in a single enthalpy term, HP , and we shall indicate when the linear approximation is possible. It should also be noted that we do not need to consider a nonlinear version of the LT model19 because the latter yields results that are indistinguishable from the linear approximation, even in the case of localized kinks.39 In order to combine the LT model with atomistic simulations, the integral in Eq. (2) must be discretized. The dislocation line is thus decomposed into elementary segments of length b, the position of which, {Yn }, is computed using the methodology explained in Sec. III B. After discretization, Eq. (2) becomes HLT ({Yn },xz ) = b
n
HP (Yn ,xz ) +
T (Yn+1 - Yn )2 , 2b2 (3)
where the sum over subscript n implicitly accounts for the periodic boundary condition in the X direction. The Hamiltonian in Eq. (3) corresponds to the 1D Frenkel-Kontorova model if the substrate potential is sinusoidal.35 The main interest of such a discretized model is that the different terms involved
144106-2
PREDICTION OF THE KINK-PAIR FORMATION . . .
PHYSICAL REVIEW B 87, 144106 (2013)
(a)
EAM 1
0
Linear enthalpy (meV/Å)
Linear enthalpy (meV/Å)
in Eq. (3) can all be computed from atomistic simulations as described in the following section.
B. Dislocation position
(b)
0
EAM 2
xz=0 MPa
xz=0 MPa
To make a link between the LT model and the atomistic simulations, the first step is to identify the dislocation position. Different methods have been proposed, for example, by adjusting the core position through comparison with solutions of the Peierls-Nabarro model41 or of anisotropic linear elasticity.30,42,43 Here, we used a simple method proposed in Ref. 38 based on the fact that the disregistery between the two Z planes above and below the dislocation glide plane (denoted Z ± ), that is, the difference in the total displacement between these two planes, increases by b each time the dislocation changes Peierls valleys. The X coordinates are thus discretized in slabs xn = nb, with n being an integer, such that in each [111] atomic column, only one atom is found in the interval In = [xn - b/2; xn + b/2]. The position in the glide plane of segment n along the dislocation line, Yn (which can be identified as a local reaction coordinate), is then defined by a 0 0 Yn = xj - xj - xj - xj , (4) b 0 0 + -
xj In Z xj In Z
-10
xz=300 MPa
-10
xz=300 MPa
-20
xz=600 MPa
-20
xz=600 MPa
(c) 10
EAM 3
(d) 10
DFT
Linear enthalpy (meV/Å)
xz=0 MPa
Linear enthalpy (meV/Å)
xz=0 MPa
0
0
-10
xz=300 MPa
-10
xz=300 MPa
-20
xz=600 MPa
-20 1
xz=600 MPa
0 0.2 0.4 0.6 0.8
0 0.2 0.4 0.6 0.8
1
Reaction coordinate Y1/a
Reaction coordinate Y1/a
where {x 0 } is the reference configuration for the straight dislocation relaxed in the initial Peierls valley. The dislocation position is scaled by a/b in order to increase by a distance between Peierls valleys a each time the dislocation changes Peierls valleys. A reaction coordinate varying between 0 and 1 can be obtained by considering Yn /a .
C. Substrate enthalpy
FIG. 1. (Color online) Minimum enthalpy path per unit length (linear enthalpy) for a straight screw dislocation computed with the NEB method and (a)(c) three different EAM potentials16,28,29 (symbols) and (d) the DFT SIESTA code. The lines are the cubic splines used to define the enthalpy function HP (Y,xz ) in Eq. (3). The splines in (c) at finite applied stress were obtained assuming that the stress dependence of the Peierls potential could be neglected.
The substrate enthalpy HP , which includes the contributions of both the Peierls potential and the work of the applied stress, was computed at the atomistic level from the hypothetical situation, where a straight infinite dislocation moves between two Peierls valleys. The substrate enthalpy HP is identified with the dislocation minimum enthalpy path (MEP) between Peierls valleys, determined here using the NEB method.44 For the three EAM interatomic potentials, the calculations were performed in a cell of width 1b, i.e., only one repeating cell in the X direction, which we call a 1b slab cell. The ° respectively. lengths in the Y and Z directions are 84 and 65 A, These dimensions are the minimum below which finite cell size effects affect significantly the dislocation energy. The number of atoms is then Nat = 900. In the 1b slab cell, translational symmetry along the X direction is imposed by means of periodic boundary conditions, such that the dislocation is constrained to glide within a straight configuration. All atoms situated in a given [111] atomic column are thus forced to move rigidly, with exactly the same displacement. The NEB implementation used here is the same as in Refs. 38 and 16, with 60 replicas and a spring constant of 1 eV/nm. The resulting enthalpies, scaled by b to obtain energies per unit length consistent with the definition in Eq. (3), are reported in Figs. 1(a)1(c) (symbols) for the three different EAM potentials as a function of Y1 , the dislocation
position along the NEB path (the subscript 1 is because only one slab in direction X is used in the simulation cell). As mentioned above, these functions include the work of the applied stress. To reconstruct the Peierls potential VP we need to correct the enthalpy for the work of the applied stress: VP = HP (y,xz ) + xz bY1 . The result is shown in Figs. 2(a)2(c) for the three EAM potentials at two different applied stresses. One can notice that the shape of the barrier is modified under the effect of the applied stress with EAM 1 and EAM 2 and not with EAM 3. More importantly, it is well known that the dislocation Peierls stress can be estimated from the maximum slope of the Peierls potential. We reported in Fig. 2 the estimated Peierls stresses for each potential and each stress (see values near triangles in Fig. 2), compared with the actual Peierls stresses P predicted by the EAM potentials. We can see that with EAM 1 and EAM 2, the estimated Peierls stress increases with the applied stress towards the true Peierls stress and reaches the Peierls stress only when the latter is applied to the system. As already observed for edge Lomer dislocations in fcc Al,38 the Peierls potential stiffens with increasing applied resolved shear stress. As a result, the Peierls stress computed from the maximum slope of the zero-stress Peierls energy will be underestimated by as much as 200 MPa in the case of potentials EAM 1 and EAM 2. Interestingly, with EAM 3, no stress dependence on the estimated Peierls stress is observed, meaning that HP is very close to the linear expansion VP (aY1 ) - xz bY1 and that the zero-stress estimate
144106-3
LAURENT PROVILLE, LISA VENTELON, AND DAVID RODNEY
(a) 12
PHYSICAL REVIEW B 87, 144106 (2013)
EAM 2
900 MPa
EAM 1 (b) 12
1250 MPa
(c) 12
1000 MPa xz=0 MPa
EAM 3
Linear energy (meV/Å)
Linear energy (meV/Å)
10 8 6 4 2 0
1090 MPa 1000 MPa
10 8 6 4 2 0
Linear energy (meV/Å)
10 8
1000 MPa
xz=0 MPa
xz=0 MPa
830 MPa 700 MPa
6
1000 MPa
4 2 0
xz=600 MPa
xz=600 MPa
xz=600 MPa
0
0.5
1
0
0.5
1
0
0.5
1
Reaction coordinate Y1/a
Reaction coordinate Y1/a
Reaction coordinate Y1/a
FIG. 2. (Color online) Peierls potentials at different applied stresses reconstructed by correcting the linear enthalpies of Fig. 1 for the work of the applied stress. The potentials were shifted in both reaction coordinate and energy in order to overlap the potential minima at each applied stress. The triangles indicate the maximum slope in each curve.
of the Peierls stress is accurate for EAM 3, in contrast to potentials EAM 1 and EAM 2. In Fig. 2(c), the only difference between the substrate potentials at different applied stresses remains due to the different equilibrium positions, and a simple translation in the glide direction suffices to superimpose the two curves. The DFT calculations were performed in a 1b slab cell containing 273 atoms within the quadrupolar arrangement. ° respectively. Lengths in the Y and Z directions are 49 and 26 A, The reaction coordinate method (or drag method)45 is used to determine the Peierls barrier with 11 images. This method was shown in Ref. 43 to be equivalent to the NEB method for the simple path considered here and is preferred for DFT calculations because of its smaller computational cost. The results are reported in Fig. 1(d). It is worth noting that the shape of the Peierls barrier obtained with EAM 3 is similar to that obtained with DFT, i.e., with no intermediate metastable state and with similar amplitudes. This is related to the detailed trajectory of the dislocation between Peierls valleys,30 which does not pass through the split configuration with EAM 3 and DFT, at variance with EAM 1 and EAM 2 (see Ref. 43 for details). On the basis of the similarities between EAM 3 and DFT, we assume that the stress dependence of the substrate enthalpy is linear in DFT, as in EAM 3. In Fig. 1(d), solid lines correspond to the results obtained from the calculation of HP through the linear approximation VP (Y1 ) - xz bY1 for different stresses.
D. Line tension
To compute the line tension T a dislocation bow-out has to be included, and thus the translational symmetry along the dislocation line must be broken. To this end we considered a 2b slab cell, i.e., a cell with a length 2b in the X direction. With EAM potentials, the minimum dimensions in the Y and ° Z directions to minimize finite-size effects are 70 and 40 A, respectively. The number of atoms is then Nat = 780. When a straight dislocation changes Peierls valleys, the three [111] columns directly in contact with the dislocation core [see Fig. 3(a)], where helicity of the bcc {111} stacking is reversed, have much larger displacements than the rest of the crystal. Atoms in these columns are identified in Fig. 3(a) by labels A
to F. Atoms A to C are in a (111) slab of width b; atoms D to F are the corresponding atoms in the second slab constituting the simulation cell. By periodicity, the crystal repeats itself, and the next atoms in the (111) stacking would be periodic images of atoms A to C. When the dislocation advances by one Peierls valley to the right, atoms A and D (above the glide plane) move in the positive X direction, while atoms B and C, and E and F (below the glide plane), move in the opposite X direction with a displacement equal to half that of atoms A and D, in absolute value [see Fig. 3(b)]. During the transition, ° the total displacement of atoms A and D is b/3 0.8 A. To compute the line tension, we have to estimate the energetic cost of bowing out the dislocation in a way consistent with the formation of a kink pair. In the latter process, a dislocation segment starts to glide towards the next Peierls valley while the rest of the dislocation remains in the initial Peierls valley. To mimic this process, we considered a dislocation relaxed under zero stress and constrained atoms A, B, and C to follow the trajectory described above: we imposed a displacement dx to atom A and -dx/2 to atoms B and C, while the atoms in the other slab, D, E, and F, were fixed at their initial positions along the X direction. We thus imposed different dislocation positions in the two slabs, Y1 0, Y2 2dx , modeling in this way the very first stage of kink-pair formation. We increased dx by increments, and between each increment, the simulation cell energy was minimized, freezing the displacements along the X direction of atoms A to F, with all other degrees of freedom being unconstrained. The energy difference E = Vat ({ri }) - Vat ({ri0 }) and the position difference Y2 - Y1 were then calculated. According to Eq. (3) and accounting for periodic boundary condition along the X direction, the quadratic term T /b[Y2 - Y1 ]2 is equal to the difference between E and the Peierls energies of the two dislocation segments. The function ELT = E - bHP (Y1 ,0) - bHP (Y2 ,0) (5)
is therefore a quadratic function of [Y2 - Y1 ] with a curvature equal to T /b, allowing us to deduce the line tension T . We report the evolution of ELT as a function of (Y2 - Y1 )/a as symbols in Fig. 3(c). The quadratic fits are also represented as solid lines. The fits are very satisfactory for the
144106-4
PREDICTION OF THE KINK-PAIR FORMATION . . .
PHYSICAL REVIEW B 87, 144106 (2013)
FIG. 3. (Color online) (a) View of the simulation cell with a 2b slab in the X direction, showing only the atoms in the dislocation core. The labels (AF) and the colors [blue (dark gray) and orange (light gray)] distinguish the atoms with larger displacements when the dislocation changes Peierls valleys (see text for details). The crystal planes Z + and Z - are also indicated. (b) Displacements of the atoms identified in (a) along the minimum energy path of a straight dislocation in the absence of applied stress for different EAM potentials. (c) Elastic energy of the dislocation computed from atomistic simulations in the 2b slab cell (symbols) and from a LT model with different adjusted values of the line ° T (EAM 2) = 3.0 eV/A, ° T (EAM 3) = 1.6 eV/A, ° and T (DFT) = 4.1 eV/A. ° tension T : T (EAM 1) = 3.8 eV/A,
three EAM potentials in the interval [Y2 - Y1 ] [0,0.05], and ° T (EAM 2) = 3.0 eV/A, ° we obtained T (EAM 1) = 3.8 eV/A, ° Beyond a position difference of and T (EAM 3) = 1.6 eV/A. 0.05, nonlinear contributions appear in ELT , which would require a higher-order expansion in Eq. (5). However, we checked that accounting for such nonlinear terms does not change significantly the kink-pair formation enthalpy, and we will consider here only the linear approximation. We performed the same type of calculations using DFT. The simulation cell was also a 2b slab in the X direction, with the quadrupolar geometry as in Sec. III C. The dimensions in ° with a total number the Y and Z directions were 37 and 18 A, of atoms of Nat = 270. The results for ELT [Eq. (5)] are shown in Fig. 3(c) as circles, together with the corresponding ° which is close quadratic fit. We obtained T (DFT) = 4.1 eV/A, to the line tension obtained with EAM 1, while, surprisingly, potential EAM 3, which has a Peierls potential very close to the DFT Peierls potential, has a line tension significantly smaller than DFT. This stresses the importance of atomic-scale effects on the dislocation line tension.
IV. CALCULATION OF THE KINK PAIR
The kink-pair enthalpy can be computed directly in large simulation cells if we employ interatomic potentials. We used a simulation cell of width 54b in the direction of the dislocation line. The dimensions in the Y and Z directions are 139 and ° respectively, for a total number of atoms of Nat = 88 A, 123 120. Again, these dimensions are the minimum to reduce finite-size effects on kink-pair formation. The dislocation MEPs are computed at different applied stresses with the NEB
method.38 Starting from initial paths that contain an expanding kink pair in order to break the translational symmetry along the dislocation line, the enthalpy variation along the MEPs is obtained from Eq. (1), and the shape of the dislocation in each NEB replica {Y (k )n } is obtained from Eq. (4). The saddle state corresponds to the replica of maximum enthalpy along the MEP, and the associated enthalpy variation corresponds to the kink-pair formation enthalpy. The results obtained with the three EAM potentials are reported in Figs. 4(a)4(c) . In these figures, we also reported the corresponding kink-pair enthalpy obtained from the LT model with a line tension fixed to the values found in Sec. III D and for a line length also equal to 54b. We should note that no fitting parameter is used here since both ingredients of the model, the line tension and the Peierls enthalpy, were taken directly from atomistic simulations. The kink-pair configurations were obtained using a Newton-Raphson method starting the computation with a trial function made of a combination of hyperbolic tangents.38 Excellent agreement is found for all three EAM models, thereby confirming that the calculation of the dislocation line tension correctly matches the elastic force involved in the kink-pair formation. With EAM 1, a transition occurs in the kink-pair formation enthalpy at xz 120 MPa. This transition (see the detailed study from Gordon et al.13 ) is a consequence of the double-hump profile of the Peierls barrier. Interestingly, this transition is also well reproduced by the LT model. This excellent agreement between LT predictions and direct atomistic calculations for all three interatomic potentials shows the validity of our approach and, in particular, that the possible variation of the line tension between the bottom and top of the Peierls valley can be neglected. We may therefore use the
144106-5
LAURENT PROVILLE, LISA VENTELON, AND DAVID RODNEY
(a)
0.8
PHYSICAL REVIEW B 87, 144106 (2013)
(b)
0.8
Atomistic simulations Line tension (EAM 1)
Kink-pair enthalpy (eV)
Kink-pair enthalpy (eV)
Atomistic simulations Line tension (EAM 2)
° Interestingly, this width depends only moderately on 40 A. the applied stress, while the height of the kink pair decreases rapidly when the applied stress increases.
0.6
0.6
EAM 1
0.4
EAM 2
0.4
V. DISCUSSION
0.2
0.2
0 0
Applied stress (MPa)
500
1000
0 0
Applied stress (MPa)
Kink-pair energy [ 0.86 eV ]
Line tension (DFT)
500
1000
(c)
0.8
(d)
0.8
Atomistic simulations Line tension (EAM 3)
0.6
0.6
EAM 3
0.4
DFT
0.4
0.2
0.2
Peierls stress [1250 MPa]
0 0
Applied stress (MPa)
500
1000
0 0
Applied stress (MPa)
500
1000
FIG. 4. (Color online) (a)(c) Kink-pair formation enthalpy computed from atomistic simulations (open symbols) compared to LT predictions (lines with solid symbols) for the same dislocation length. (d) LT model prediction based on DFT data.
DFT line tension and Peierls barrier to predict the kink-pair formation enthalpy from DFT calculations. The results are reported in Fig. 4(d) and are discussed in the following section. From these calculations, we can also obtain the profile of the critical kink pair on the screw dislocation at different applied stresses (Fig. 5). The width of the critical kink pair is about
1.5
Kink-pair profile (Å)
100 MPa 300 MPa 600 MPa
1
0.5
0 -60 -40 -20 0 20 40 60 80 100
Dislocation line (Å)
FIG. 5. (Color online) Profiles of the dislocation with a critical kink pair at different applied stresses, computed from the LT model with parameters from DFT calculations.
The dislocation line tension was calculated from atomicscale simulations, with the in-row atomic displacements constrained to reproduce the first stage of the formation of a kink pair on the dislocation. The present work is readily applicable to screw dislocations in all bcc metals and paves the way for extensions to other dislocations in more complex materials. Once parameterized, the LT model allows us to predict very accurately the kink-pair formation enthalpy computed with three different EAM potentials that have distinct Peierls barriers and line tensions. On the basis of this validation, we employed the same method to predict the kink-pair formation enthalpy from DFT calculations. This yields a DFT Peierls stress of about 1.25 GPa, as reported in Fig. 4(d), in agreement with previous direct calculations using DFT.11,43 The present Peierls stress computations ignored the effect of zero-point vibrations that must be accounted for to compare with experimental data.16 The kink-pair formation energy, i.e., the kink-pair formation enthalpy at zero applied stress, is found to be Ekp = 0.86 eV, close to the experimental estimates, which lie between Ekp = 0.6 and 0.91 eV, depending on the temperature regime considered.46,47 This calculation will have to be confirmed using different DFT codes, as was performed for the screw dislocation Peierls barrier in Ref. 43. With EAM, the kink-pair formation energy varies from Ekp = 0.7 eV for EAM 1 to Ekp = 0.55 eV for EAM 3. In direct computations of kink-pair formation with EAM potentials, we have verified that the formation energy is not impacted by an increase of the simulation cell in the X direction, i.e., the dislocation length, meaning that the interaction between the kink pairs due to periodic boundary conditions is negligible in the present geometry. Recently, the kink-pair formation energy was estimated using DFT by Itakura et al.,30 who found Ekp = 0.73 eV, using a different method of calculation and a different DFT code (VASP). These authors also estimated the ° screw dislocation line tension and found T (EAM) = 3.8 eV/A using EAM variant 5 proposed by Mendelev et al.,28 while ° using the DFT VASP code. These results T (DFT) = 3.08 eV/A are consistent with the results presented here. In the method proposed by Itakura et al., only one atom (not specified by these authors) was displaced, while according to our calculations, a satisfactory precision on the line tension requires us to constrain the displacement of all A, B, and C atoms (see Fig. 3). Moreover, our study shows that the line tension computation needs only a slab of width 2b, while in the above work, a slab of width 3b was used. This allows us to reduce the required number of atoms, of marked interest for DFT calculations. As mentioned in the Introduction, the linear elasticity calculation is based on the hypothesis that the distance over which the dislocation bows out is large compared to the amplitude of the bow-out, an assumption which is not adapted to kink pairs. In the literature, the line tension is often estimated
Kink-pair enthalpy (eV)
Kink-pair enthalpy (eV)
144106-6
PREDICTION OF THE KINK-PAIR FORMATION . . .
PHYSICAL REVIEW B 87, 144106 (2013)
from the elasticity theory of dislocations18,23,24 and contains both elastic and core terms. The elastic term is expressed for an isotropic linear elastic medium as Tel = b2 /4 (1 - ) ln(R/r0 )[(1 + ) cos2 + (1 - 2 ) sin2 ], where is the shear modulus, is Poisson's ratio, b is the dislocation Burgers vector, is the angle between b and the dislocation line direction, and r0 and R are the inner and outer cutoff radii, respectively. When treating specific dislocation configurations, the logarithmic term systematically depends on the length of the bowing segment instead of the outer cutoff radius R .24 The above expression is thus an average estimate, which assumes that the dislocation bows out between obstacles made of dislocation junctions, such that if is the dislocation density, we have the usual estimate R 1/ . This assumption is well adapted to fcc crystals where the Peierls valleys are negligible and the dislocations smoothly bow out between obstacles. In this case, the line tension depends on the distance between dislocations and, in the context of atomistic simulations, depends logarithmically on the simulation cell size.5,38 On the other hand, in the case of kinks, the bowing segments are very short (see Fig. 5) and are not related to the dislocation density or the simulation cell size. In the case of kinks, we therefore expect that the line tension depends on the length of the bowing segment , rather than the simulation cell size. This is indeed what we observed with the EAM potentials. Varying the cell size by one order of magnitude, we did not find any variation of our estimate of the line tension T . A consequence is that the LT prediction for the kink-pair formation enthalpy is independent of the cell size, a property that we also checked by direct NEB calculations in cells of varying cell dimensions. Also, we should note that in our calculations, the elastic and core terms cannot be computed separately, and since r0 , an estimate of the elastic part based on the elasticity theory would be highly uncertain. The only conclusion we can draw is that, since all potentials predict the same elastic constants, if elasticity dominated the line tension, then all potentials should have the same line tension. Since we found a large variability between potentials, we conclude that the core contribution to the line tension is large, at least for kinks on bcc screw dislocations. For researchers who employ the above formula for Tel , it is of some interest to provide a quantitative estimate. Using the standard elastic coefficients for iron, i.e., = (C11 - C12 + C44 )/3 = 71.5 GPa and = 0.29 and 40 ° i.e., the half width of a kink pair (see Fig. 5), and r0 b, A,
° which is by far smaller than values we found Tel = 1.05 eV/A, worked out from atomic-scale simulations. One last point of interest in the comparison between LT model and atomic-scale simulations has a bearing on the dynamics of the dislocation, which is briefly discussed here. The time derivation of reaction coordinates {Yn } in Eq. (3), in which kinetic energy can be easily added, leads to a dynamical equation with soliton solutions.48 Such solitons actually correspond to traveling kinks. In the LT model, the peculiarity of solitons is that their collision can be conservative, i.e., when a left kink encounters a right kink, a new kink pair is created spontaneously where the two kinks met, allowing the elastic line to cross the next Peierls hill. Said in other words, within the LT model, two kinks can cross each other without annihilation, leading to avalanches in the dislocation dynamics where several Peierls valleys can be overcome in a single dynamical event. When the kink number is small, their collisions are solely due to the periodic boundary conditions, i.e., it is a computational artifact if one considers how plastic deformation truly proceeds.49 To determine whether such soliton solutions exist in atomic-scale simulations,7 the trajectory of a screw dislocation was simulated in the 54b cell, with molecular dynamics (MD) in the NVE microcanonical statistical ensemble, i.e., with constant particle number, volume and energy. The starting state is the kink-pair saddle state computed from NEB with velocities randomly distributed with a zero-centered Gaussian distribution normalized to fulfill the condition on total energy. The same type of simulation was repeated for different stresses, ranging from 200 to 800 MPa, and for different total energies. In spite of many attempts, no soliton-like solution emerged from the collision between kinks. They annihilate systematically when colliding after traveling along the dislocation segments in opposite directions and passing through the periodic boundary conditions. Thus it seems that the interaction of the traveling kinks with the crystal vibrational modes36,50 is sufficiently important in a real three-dimensional crystal to cancel the solitary waves, in contrast to the simplified two-dimensional space considered in the LT model, a point which deserves proper study at the atomistic level.
ACKNOWLEDGMENTS
The authors thank Franc ¸ ois Willaime for fruitful discussions.
1 2
L. P. Kubin, Rev. Deform. Behav. Mater. 1, 244 (1977). D. Caillard and J. L. Martin, Thermally Activated Mechanisms of Crystal Plasticity (Pergamon, Amsterdam, 2003). 3 S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett. 84, 1499 (2000). 4 A. H. W. Ngan and M. Wen, Phys. Rev. Lett. 87, 075505 (2001). 5 L. H. Yang, P. Soderlind, and J. A. Moriarty, Philos. Mag. A 81, 1355 (2001). 6 S. L. Frederiksen and K. W. Jacobsen, Philos. Mag. 83, 365 (2003). 7 J. Chaussidon, M. Fivel, and D. Rodney, Acta Mater. 54, 3407 (2006).
L. Ventelon, F. Willaime, and P. Leyronnas, J. Nucl. Mater. 386-388, 26 (2009). 9 C. Woodward and S. I. Rao, Phys. Rev. Lett. 88, 216402 (2002). 10 L. Ventelon and F. Willaime, J. Comput. Aided Mater. Des. 14, 85 (2007). 11 L. Ventelon and F. Willaime, Philos. Mag. 90, 1063 (2010). 12 E. Clouet, L. Ventelon, and F. Willaime, Phys. Rev. B 84, 224107 (2011). 13 P. A. Gordon, T. Neeraj, Y. Li, and J. Li, Modell. Simul. Mater. Sci. Eng. 18, 085008 (2010). 14 C. Domain and G. Monnet, Phys. Rev. Lett. 95, 215506 (2005).
8
144106-7
LAURENT PROVILLE, LISA VENTELON, AND DAVID RODNEY
15
PHYSICAL REVIEW B 87, 144106 (2013)
32 33
L. Romaner, C. Ambrosch-Draxl, and R. Pippan, Phys. Rev. Lett. 104, 195503 (2010). 16 L. Proville, D. Rodney, and M.-C. Marinica, Nat. Mater. 11, 845 (2012). 17 E. Orowan, Proc. Phys. Soc. 52, 8 (1940). 18 A. Seeger and U. Holzwarth, Philos. Mag. 86, 3861 (2006). 19 J. E. Dorn and S. Rajnak, Trans. AIME 230, 1052 (1964). 20 J. Friedel, Dislocations (Pergamon, London, 1964). 21 A. Seeger and P. Schiller, in Physical Acoustics, edited by W. Mason (Academic, New York, 1966), Vol. 3A, p. 361. 22 A. Seeger, in Dislocation 1984, edited by P. Veyssi` ere, L. Kubin, and J. Castaing (CNRS, Paris, 1984). 23 J. P. Hirth and J. Lothe, Theory of Dislocations (Wiley, New York, 1982). 24 J. P. J. P. Hirth, T. Jossang, and J. Lothe, J. Appl. Phys. 37, 110 (1966). 25 J. M. Soler, E. Artacho, J. D. Gale, A. Garc´ ia, J. Junquera, P. Ordej´ on, and D. S´ anchez-Portal, J. Phys. Condens. Matter. 14, 2745 (2002). 26 P. Guyot and J. E. Dorn, Can. J. Phys. 45, 983 (1967). 27 S. B. Sinnott and S. W. Brenner, MRS Bull. 37, 469 (2012). 28 M. I. Mendelev, S. W. Han, D. J. Srolovitz, G. J. Ackland, D. Y. Sun, and M. Asta, Philos. Mag. 83, 3977 (2003). 29 P. A. Gordon, T. Neeraj, and M. I. Mendelev, Philos. Mag. 91, 3931 (2011). 30 M. Itakura, H. Kaburaki, and M. Yamaguchi, Acta Mater. 60, 3698 (2012). 31 D. J. Bacon, Y. N. Osetsky, and D. Rodney, in Dislocations in Solids, edited by L. K. J. P. Hirth (Elsevier, Amsterdam, 2008), p. 1.
D. Rodney, Phys. Rev. B 76, 144108 (2007). J. Frenkel and T. Kontorova, Phys. Z. Sowj. 13, 1 (1938) [J. Phys. (Moscow) 1, 137 (1939)]. 34 B. Jo´ os, Solid State Comm. 42, 709 (1982). 35 B. Jo´ os and M. S. Duesbery, Phys. Rev. B 55, 11161 (1997). 36 J. D. Eshelby, Proc. R. Soc. London, Ser. A 266, 222 (1962). 37 K. Kanga, V. V. Bulatov, and W. Cai, Proc. Natl. Acad. Sci. USA 109, 15174 (2012). 38 D. Rodney and L. Proville, Phys. Rev. B 79, 094108 (2009) 39 D. Rodney and J. Bonneville, in Physical Metallurgy, edited by K. Hono and D. Laughlin (Elsevier, Amsterdam, 2013). 40 R. Gr¨ oger and V. Vitek, Acta Mater. 56, 5426 (2008). 41 L. Pizzagalli, A. Pedersen, A. Arnaldsson, H. J´ onsson, and P. Beauchamp, Phys. Rev. B 77, 064106 (2008). 42 R. Gr¨ oger and V. Vitek, Modell. Simul. Mater. Sci. Eng. 20, 035019 (2012). 43 L. Ventelon, E. Clouet, F. Willaime, and D. Rodney, Acta Mater. (in press). 44 H. J´ onsson, G. Mills, and K. Jacobsen, in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by G. C. B. J. Berne and D. F. Coker (World Scientific, Singapore, 1998), p. 385. 45 G. Henkelman, G. J´ ohannesson, and H. J´ onsson, in Progress on Theoretical Chemistry and Physics (Kluwer, Dordrecht, 2000), p. 269. 46 D. Brunner and J. Diehl, Phys. Status Solidi A 124, 155 (1991). 47 D. Brunner and J. Diehl, Phys. Status Solidi A 124, 455 (1991). 48 M. Peyrard and T. Dauxois, Physics of Solitons (Cambridge University Press, Cambridge, 2006). 49 D. Caillard, Acta Mater. 58, 3504 (2010). 50 A. D. Brailsford, J. Appl. Phys. 41, 4439 (1970).
144106-8