Screw dislocation in zirconium: An ab initio study
aa r X i v : . [ c ond - m a t . m t r l - s c i ] O c t Screw dislocation in zirconium: an ab initio study
Emmanuel Clouet ∗ CEA, DEN, Service de Recherches de M´etallurgie Physique, 91191 Gif-sur-Yvette, France (Dated: November 7, 2018)Plasticity in zirconium is controlled by 1 / h i screw dislocations gliding in the prism planesof the hexagonal close-packed structure. This prismatic and not basal glide is observed for a givenset of transition metals like zirconium and is known to be related to the number of valence electronsin the d band. We use ab initio calculations based on the density functional theory to study thecore structure of screw dislocations in zirconium. Dislocations are found to dissociate in the prismplane in two partial dislocations, each with a pure screw character. Ab initio calculations also showthat the dissociation in the basal plane is unstable. We calculate then the Peierls barrier for ascrew dislocation gliding in the prism plane and obtain a small barrier. The Peierls stress deducedfrom this barrier is lower than 21 MPa, which is in agreement with experimental data. The abilityof an empirical potential relying on the embedded atom method (EAM) to model dislocations inzirconium is also tested against these ab initio calculations.
PACS numbers: 61.72.Lk, 61.72.Bb
I. INTRODUCTION
Plasticity in α -zirconium is controlled by dislocationswith a 1 / h i Burgers vector gliding in the prismplanes of the hexagonal compact (hcp) lattice.
The rel-ative ease of prismatic glide compared to basal glide hasbeen shown to be linked to the ratio of the correspond-ing stacking fault energies, which in turn is controlled bythe electronic structure. In particular, Legrand used atight binding model to show that prismatic slip in tran-sition metals of the IV B column (Zr, Ti, Hf) originatesfrom the electronic filling of the valence d band. As aconsequence, it appears necessary to take into accountthe anisotropy of the d orbital, and hence the angulardependence of the atomic bonding, to model dislocationsin these transition metals . One cannot therefore rely oncentral forces empirical potential and needs to take ac-count of the electronic structure. Tight binding models or ab initio calculations show indeed that a 1 / h i screw dislocation, either in Zr or in Ti, spreads in prismplanes, in agreement with the prismatic glide observedexperimentally. But none of these atomistic simulationscalculate the Peierls stress of a screw dislocation. Ac-cording to some authors, its core structure is not com-pletely planar, which may be the cause of a high Peierlsstress.It is true, experimentally, that screw dislocationsglide with difficulty compared to other dislocation char-acters in zirconium or titanium alloys: characteris-tic microstructures, with long and straight dislocationsaligned along their screw orientation, are observed atlow temperature, and the flow stress is stronglytemperature dependent, in agreement with theassumption of a high Peierls barrier which must be over-come by the nucleation of double kinks. But experimentsalso show that the yield stress in zirconium or in tita-nium strongly decreases with a decreasing amount of in-terstitial impurities like oxygen. It is there-fore probable that the high Peierls stress of screw dis- locations has an extrinsic cause. In pure zirconium orpure titanium, this Peierls stress may not be as high andscrew dislocations are probably gliding as easily as otherorientations.Recently, Mendelev and Ackland developed an em-pirical potential for Zr in the Embedded Atom Method(EAM) formalism. Using this potential, Khater andBacon showed that it leads to a screw dislocation thatspontaneously spreads in the prism plane, the configu-ration dissociated in the basal plane being metastable.They also showed that the Peierls stress of a screw dislo-cation gliding in the prism plane is not so different fromthe Peierls stress of an edge dislocation and that thisstress is small (22 MPa for the screw and 16 MPa for theedge). These results therefore support experimental find-ings stating that screw dislocations are gliding in prismplanes with a low Peierls stress in pure zirconium. Butthese simulations rely on a central forces potential, whichis not well suited to describe dislocations in a hcp tran-sition metal like Zr, as described above. More reliableatomistic simulations, incorporating a description of theelectronic structure, are therefore needed to confirm thiseasy glide of screw dislocations in pure zirconium.This article aims to use ab initio calculations so asto fully characterize the core structure of a 1 / h i screw dislocation in zirconium and estimate its Peierlsstress. We also examine generalized stacking faults asdislocation core structures are closely related to them.Two different ab initio methods, Siesta and Pwscf ,have been used. Siesta offers the advantage of effi-ciency, allowing simulating more atoms than with a stan-dard ab initio code, whereas
Pwscf offers the advantageof robustness. All calculations are also performed withMendelev and Ackland EAM potential, so as to iden-tify its ability to model dislocations in zirconium. Inaddition, this empirical potential is used to study theconvergence of our results with the size of the simulationcell. II. ATOMIC INTERACTION MODELING
Atomistic simulations have been performed both withan empirical interatomic potential and with ab initio cal-culations. The empirical potential that we used is theEAM potential developed by Mendelev and Ackland. This potential is labeled ab initio values of the stacking fault energies in the basal and prismplanes have been included in the fitting procedure. Us-ing this potential, Khater and Bacon showed that a1 / h i screw dislocation spontaneously dissociates inthe prism plane and that a metastable configuration dis-sociated in the basal plane also exists.The ab initio calculations are relying on the DensityFunctional Theory (DFT) in the Generalized GradientApproximation (GGA) with the functional proposed byPerdew, Burke, and Ernzerhof (PBE) and the pseudopo-tential approximation. Two different ab initio codes areused, Siesta and Pwscf . In the
Siesta code, valence electrons are describedby a localized basis set corresponding to a linear com-bination of pseudoatomic orbitals with 13 functions peratom. We used a norm conserving pseudopotential ofTroulliers-Martins type with 4p electrons included assemicore. Electronic density of state is broaden with theMethfessel-Paxton function with a broadening of 0.3 eVand the integration is performed on a regular grid of14 × × i.e. the basis and the pseu-dopotential, has already been used to study vacancy dif-fusion in zirconium and comparison with plane wavesDFT calculations has led to a reasonable agreement.In the Pwscf code, valence electrons are describedwith plane waves using a cutoff energy of 28 Ry. Thepseudopotential is ultrasoft of Vanderbilt type with 4sand 4p electrons included as semicore. The same k-point grid and the same electronic broadening are usedas with
Siesta code.To validate these different atomic interaction models,it is worth comparing their results to available experi-mental data for some bulk properties of Zr. All modelslead to an equilibrium lattice parameter and a c/a ratioin good agreement with experimental data (Tab. I). Inparticular, a ratio lower than the ideal p / ∼ . as the hcp lattice contains two atoms in its primi-tive unit cell, some internal degrees of freedom may existwhen applying a homogeneous strain. One needs to allowfor atomic relaxations when computing C , C or C constants. It is also possible to calculate inner elastic
TABLE I. Bulk properties of hcp Zr calculated with differ-ent atomic interaction models and compared to experimentaldata: lattice parameter a , c/a ratio of the hexagonal lattice,relaxed elastic constants C ij , inner elastic constants e ij and d ij , phonon frequencies ω and ω of the optical branches atthe Γ point (Eq. 1), inner elasticity contribution to elasticconstant δC (Eq. 2).Expt. EAM Siesta Pwscf a (˚A) 3.232 c/a C (GPa) 155.4 a C (GPa) 172.5 a C (GPa) 67.2 a
75. 86. 70. C (GPa) 64.6 a
76. 68. 65. C (GPa) 36.3 a
44. 24. 26. C (GPa) 44.1 a e (meV ˚A − ) 27.0 17.0 18.6 e (meV ˚A − ) 118. 122. 101. d (meV ˚A − ) 30.0 38.9 36.6 ω (THz) 2 . ± . ω (THz) 4 . ± . δC (GPa) 5.33 14.3 11.5 a Experimental elastic constants have been measured at 4 K. constants to characterize these internal degrees of free-dom. These are also given in Tab. I using the notationsintroduced by Cousins. Two of these inner elastic con-stants, e and e , are related to the phonon frequenciesof the optical branches at the Γ point ω i = 2 r Ω e ii m , (1)where Ω = a c √ / m theatomic mass. The last inner elastic constants d cou-ples the internal degrees of freedom to the homogeneousstrain. It leads to a contribution to the elastic constants δC = d e . (2)If C ij are the unrelaxed elastic constants, i.e. the elasticconstants calculated by imposing a homogeneous strainto the hcp lattice without letting atoms relax from theirinitial positions, the true elastic constants are given by C = C − δC , C = C + δC , and C = C − δC , all other elastic constants being unchanged. III. STACKING FAULT ENERGIES
Dislocation dissociation is controlled by the existenceof a metastable stacking fault for the correspondingplane. According to the results obtained by Khater andBacon with Mendelev and Ackland EAM potential, TABLE II. Stacking fault energies in the basal plane, γ b , andin the prism plane, γ p , calculated with different atomic inter-action models, including Vasp calculations of Domain et al .. and of Udagawa et al .. R is the ratio defined by Legrand, and d eqb and d eqp are the dissociation lengths of a screw dislo-cation respectively in the basal (Eq. 5) and in the prism plane(Eq. 6). EAM Siesta Pwscf Vasp Vasp γ b (mJ m − ) 198. 199. 213. 200. 227. γ p (mJ m − ) 135. (274.) a R = C γ b /C γ p d eqb (˚A) 4.0 2.0 2.7 3.4 3.2 d eqp (˚A) 7.8 4.6 5.9 9.6 7.4 a For the EAM calculation of the prismatic stacking fault energy,the value in parenthesis corresponds to the maximum in a/ a/ . c [0001]. a 1 / h i screw dislocation can dissociate either in abasal or in a prism plane. To characterize these even-tual dissociations, we compute generalized stacking faultenergies – or γ -surfaces – for both the basal and theprism planes. A. Methodology γ -surfaces describe the energy variation when twoparts of a crystal are rigidly shifted for different faultvectors lying in a given crystallographic plane. Atomsare allowed to relax in the direction perpendicular to thefault plane. We calculate these γ -surfaces for both thebasal and prism planes using full periodic boundary con-ditions. To introduce only one fault in the simulationcell, the same shift as the one corresponding to the faultvector is applied to the periodicity vector perpendicularto the fault plane. No free surfaces is therefore intro-duced in the simulation cell, which allows a fast conver-gence of the fault energies with the number of stackedplanes. A periodic stacking of at least 16 { } planesis used for the basal fault and 12 { } planes for theprismatic fault. This corresponds to a distance betweenfault planes h = 8 c ∼
41 ˚A and h = 6 √ a ∼
34 ˚A.Increasing the number of planes in the stacking modifiesthe energies by less than 1 mJ m − . Generalized stackingfault energies are calculated on a regular grid of 10 × FIG. 1. Generalized stacking fault energy in the basal planecalculated with (a)
Pwscf , (b)
Siesta and (c) EAM poten-tial. The arrows indicate Burgers vectors of the partial dis-locations corresponding to a dissociation in the basal plane.The dashed line is the [1¯100] direction used in Fig. 2. Contourlines are drawn at the base every 50 mJ m − . B. γ -surfaces
1. Basal plane γ -surfaces corresponding to the basal plane are shownin Fig. 1 for the three interaction models we used. In allcases, a local minimum can be found at 1 / I fault. This minimumdoes not vary when full atomic relaxations are allowed g ( J m − ) g ( m e V Å − ) a [11−00]PwscfSiestaEAM FIG. 2. Generalized staking fault energy in the basal planealong the [1¯100] direction ( cf . Fig. 1) calculated with Pwscf , Siesta and EAM potential. instead of being constrained to the direction perpendicu-lar to the fault plane. All methods give a close value forthe fault energy γ b in this minimum (Tab. II). A goodagreement is also obtained with the values calculated byDomain et al . and by Udagawa et al . using Vasp abinitio code. The depth of this local minimum is morepronounced with the empirical EAM potential [Fig. 1(c)]than in ab initio calculations [Fig. 1(a) and (b)]. This ap-pears clearly in Fig. 2 where the fault energy predicted bythe different interaction models are compared along the[1¯100] direction. We will see latter that this has conse-quences on the stability of a screw dislocation dissociatedin the basal plane.The γ -surface calculated with the EAM empirical po-tential has another minimum located in 2 / Ab initio calculations confirm that this faultvector gives a maximum. Finally, it is worth pointingthat both ab initio techniques give a very similar γ -surface: the shapes are identical and the amplitudes donot differ by more than 10%.
2. Prism plane γ -surfaces for the prism plane are shown in Fig. 3.Both ab initio methods lead to a similar γ -surface[Figs. 3(a) and (b)]. In particular, both Pwscf and
Siesta predicts the existence of a minimum at halfway ofthe Burgers vector, i.e. in 1 / γ -surface along the [1¯210] direction (Fig. 4), thisminimum is a little more pronounced with Pwscf thanwith
Siesta . The same minimum was also present in the
Vasp calculations of Domain et al ., but they obtained alower value γ p of the fault energy in this point (Tab. II). FIG. 3. Generalized stacking fault energy in the prism planecalculated with (a)
Pwscf , (b)
Siesta and (c) EAM poten-tial. The arrows indicate Burgers vectors of the partial dis-locations corresponding to a dissociation in the prism plane.Contour lines are drawn at the base every 75 mJ m − . On the other hand, Udagawa et al . obtained a valueclose to our result using also Vasp ab initio code withthe PAW method and the PBE-GGA functional. Theypointed out that the discrepancy arises from an insuf-ficient number of stacked planes in the γ -surface calcu-lation of Domain et al .. The energy of the metastablestacking fault energy in the prism plane appears there-fore higher than the value 145 mJ m − initially suggestedby Domain et al .: both our Pwscf and Siesta cal-culations, as well as the Udagawa et al . result, leads to g ( J m − ) g ( m e V Å − ) a/3 [12−10]PwscfSiestaEAM FIG. 4. Generalized staking fault energy in the prism planealong the [1¯210] direction calculated with
Pwscf , Siesta andEAM potential. a value of about 200 mJ m − .The γ -surface calculated with the EAM potentialis quite different [Fig. 3(c)] as the point located in1 / ab initio calculations. A minimum is found for apoint located in a/ . c [0001]. One thereforeexpects that a 1 / h i { } dislocation dissociatesinto two partial dislocations with a Burgers vector com-ponent orthogonal to the one of the perfect dislocation.In particular, a screw dislocation should dissociate in twopartial dislocations with a small edge character. Khaterand Bacon showed that the empirical potential of Ack-land et al . suffers from the same artefact. As noted byBacon and Vitek, all central forces potentials stabilizeindeed such a stacking fault a/ αc [0001] with α = 0, a minimum also predicted by a simple hard spheremodel. This minimum either disappears or is locatedexactly in a/ α = 0) only when the angular de-pendence of the atomic interactions is taken into account.The value γ p of the stacking fault energy obtained withMendelev and Ackland EAM potential for this mini-mum is much lower than our ab initio value (Tab. II).This is quite normal as Mendelev and Ackland used the ab initio value of Domain et al . to fit their potential. C. Dislocation dissociation
Before using atomistic simulations to obtain disloca-tion core structures and their associated Peierls stress, itis worth looking what can be learned from these γ -surfacecalculations. Legrand proposed a criterion based on theelastic constants and the metastable stacking fault en-ergies to determine if glide occurs in the base or prismplane in a hcp crystal. According to this criterion, pris-matic glide is favored if the ratio R = C γ b /C γ p islarger than 1. Tab II shows that this is the case for theEAM potential and the Pwscf calculation, as well as forthe
Vasp calculation of Domain et al . and of Udagawa et al .. On the other hand,
Siesta leads to a value tooclose to 1 to be able to decide between basal and pris-matic glide.One can also use dislocation elasticity theory to com-pute the dissociation distance of a dislocation both in thebasal and prism planes. According to elasticity theory,the energy variation caused by a dissociation of length d is ∆ E diss ( d ) = − b (1) i K ij b (2) j ln (cid:18) dr c (cid:19) + γd, (3)where b (1) and b (2) are the Burgers vectors of each par-tial dislocation, γ the corresponding stacking fault en-ergy, K the Stroh matrix controlling dislocation elas-tic energy, and r c the core radius. Minimization of thisenergy leads to the equilibrium dissociation length d eq = b (1) i K ij b (2) j γ . (4)When the hcp crystal is oriented with the x , y , and z axis respectively along the [10¯10], [0001], and [1¯210] di-rections, for a dislocation lying along the z direction, the K matrix is diagonal with its components given by K = 12 π (cid:0) ¯ C + C (cid:1) s C (cid:0) ¯ C − C (cid:1) C (cid:0) ¯ C + C + 2 C (cid:1) ,K = r C C K ,K = 12 π r C ( C − C ) , where ¯ C = √ C C .The γ -surface of the basal plane (Fig. 1) indicates apossible dissociation 1 / → / / / d eqb = (3 K − K ) a γ b . (5)According to the minimum of the γ -surface in theprism plane (Fig. 3), a 1 / / ± αc/a [0001]. The parameter α con-trols the position of the stacking fault minimum alongthe [0001] direction, i.e. α = 0 for Pwscf and
Siesta ,and α = 0 .
14 for the EAM potential. The dissociationlength of a screw dislocation in the prism plane is then d eqp = (cid:0) K a − α K c (cid:1) γ p . (6)The dissociation lengths d eqb and d eqp calculated fromthe elastic constants and the stacking fault energies aregiven in table II. Elastic constants predicted by theatomic interaction models are used in each case. For allenergy models, one expects a larger dissociation in theprism than in the basal plane. We will compare in thefollowing section these dissociation lengths predicted byelasticity theory with the ones observed in our atomisticsimulations of the dislocation core structure. IV. SCREW DISLOCATION COREA. Methodology
FIG. 5. Screw dislocation periodic arrangements used foratomistic simulations. ~U and ~U are the periodicity vectorsof the arrangement, and ~A the cut vector of the dislocationdipole. The thin vertical line corresponds to the prism glideplane. Our atomistic simulations of the core structure ofa screw dislocation are based on full periodic bound-ary conditions. This requires introducing a dislocationdipole in the simulation cell. Two periodic arrangementsof the dislocations have been used (Fig. 5). In the Oarrangement, dislocations with opposite Burgers vectorsare located on the same prism and basal planes, i.e. thetwo foreseen glide planes. On the other hand, only dislo-cation with the same Burgers vectors can be found on agiven prism or basal plane in the S arrangement.Periodicity vectors of the O arrangement are, beforeintroducing the dislocations, ~U = na [10¯10] − mc [0001], ~U = na [10¯10] + mc [0001], and ~U = a [1¯210]. Theintegers n and m are taken equal to keep an aspect ratioclose to a square. This simulation cell has been used for ab initio calculations with n = 4 (128 atoms), n = 5 (160atoms), and n = 6 (288 atoms).For the S arrangement, ~U = na [10¯10], ~U = mc [0001], and ~U = a [1¯210]. Ab initio calculationshave been performed with n = m and varying between n = 5 (100 atoms) and n = 8 (256 atoms).Both dislocation arrays are quadrupolar: the vec-tor joining the two dislocations composing the primitivedipole is ~D = 1 / ~U + ~U ). Because of the centrosym-metry of this arrangement and the symmetry propertiesof the Volterra elastic field, this ensures that the stresscreated by other dislocations is minimal at each disloca-tion position. The cut vector ~A defining the dislocation dipole is obtained by a π/ ~D .The dislocation dipole is introduced in the simulationcells by applying to all atoms the elastic displacementpredicted by anisotropic elasticity theory taking fullaccount of the periodic boundary conditions. A homo-geneous strain is also applied to the simulation cell so asto cancel the plastic strain introduced by the dislocationdipole and minimize the elastic energy. This strain isgiven by ε ij = − b i A j + b j A i S , where S = | ( ~U ∧ ~U ) · ~e z | is the surface of the simulationcell perpendicular to the dislocation lines. This homoge-neous strain adds some tilt components to the periodicityvectors. Atoms are then relaxed until all components ofthe atomic forces are smaller than 5 meV ˚A − for Siesta ,2 meV ˚A − for Pwscf , and 0.1 meV ˚A − for the EAMpotential. B. Core structure (a)
Pwscf[101-0][0001][12-10] (b)
Siesta (c)
EAM
FIG. 6. Differential displacement maps around one of the two1 / n = m = 7 (196 atoms). Atomsare sketched by circles with a color depending on the (1¯210)plane to which they belong. The arrow between two atomiccolumns is proportional to the [1¯210] component of the differ-ential displacement between the two atoms. Displacementssmaller than 0 . b are not shown. Crosses × correspond tothe positions of the two partial dislocations, and + to theirmiddle, i.e. the position of the total dislocation. Starting from perfect dislocations, atom relaxationleads to dislocations spread in the prism plane, whateverthe interaction model (EAM,
Siesta , or
Pwscf ) andwhatever the simulation cell used. This can be clearlyseen by plotting differential displacement maps as intro-duced by Vitek. These maps (Fig. 6) show that thestrain created by the screw dislocation spreads out inthe (10¯10) prism plane and that displacements outsidethis plane are much smaller.To characterize this spreading in the (10¯10) prismplane, we extract from our atomistic simulations the dis-registry D ( x ) created by the dislocation. This is de- d D / d x Position x in prismatic plane (Å) 0 0.5 1 1.5 D ( x ) / b PwscfSiestaEAM
FIG. 7. Disregistry D ( x ) created by the screw dislocationin the prism plane and corresponding dislocation density ρ ( x ) = ∂D ( x ) /∂x for the S periodic arrangement with 196atoms ( n = m = 7). Symbols correspond to atomistic sim-ulations and lines to the fit of the Peierls-Nabarro model tothese data. For clarity, disregistries D ( x ) have been shiftedby 0.25 between the different interaction models. fined as the displacement difference between the atomsin the plane just above and those just below the dis-location glide plane. The derivative of this function, ρ ( x ) = ∂D/∂x corresponds to the dislocation density.Fig. 7 shows the disregistry obtained for the three inter-action models. In all cases, the b discontinuity createdby the screw dislocation does not show a sharp interface,but spreads on a distance ∼
10 ˚A.Peierls and Nabarro built a model that leads toa simple expression of the disregistry. The analyticalexpression they obtained can be extended to a dissoci-ated dislocation. As suggested by the prismatic γ -surface(Fig. 3), we assume that the screw dislocation dissociatesin two equivalent partial dislocation separated by a dis-tance d . Based on the Peierls-Nabarro model, we writethe disregistry created by a single dislocation D dislo ( x ) = b π (cid:26) arctan (cid:20) x − x − d/ ζ (cid:21) + arctan (cid:20) x − x + d/ ζ (cid:21) + π/ (cid:27) , (7)where x is the dislocation position, d its dissociationlength, and ζ the spreading of each partial dislocation.We need then to take into account that we do not haveonly one dislocation on a given prism plane but a periodicarray (Fig. 5). The disregistry created by an array of period L is D L ( x ) = ∞ X n = −∞ D dislo ( x − nL )= b π arctan tan (cid:0) πL [ x − x − d/ (cid:1) tanh (cid:16) πζL (cid:17) + π (cid:22) x − x − d/ L + 12 (cid:23) + arctan tan (cid:0) πL [ x − x + d/ (cid:1) tanh (cid:16) πζL (cid:17) + π (cid:22) x − x + d/ L + 12 (cid:23)(cid:27) , (8)where ⌊·⌋ is the floor function. For the O arrangement(Fig. 5a), the disregistry in the prism plane should begiven by D ( x ) = D L ( x ) − D L ( x − L/
2) with L = 2 mc ,whereas it should be D ( x ) = D L ( x ) with L = mc for theS arrangement (Fig. 5b).We fit this analytical expression of the dislocation dis-registry to the data coming from the atomistic simula-tions. Fig. 7 shows a good agreement between atomisticsimulations and the model, using only three fitting pa-rameters: the dislocation position x , the dissociationlength d and the spreading ζ . This procedure thereforeallows us to determine the location of the dislocation cen-ter. For all interaction models, we find that this centerlies in between two (0001) atomic planes. One can see inFig. 6 that this position corresponds to a local symmetryaxis of the differential displacement map. This is differ-ent from the result obtained by Ghazisaeidi and Trinkle in Ti where the center of the screw dislocation was foundto lie exactly in one (0001) atomic plane, a position thatcorresponds in Zr to the saddle point between two Peierlsvalleys, as it will be shown below.The dissociation length of the screw dislocation ob-tained through this fitting procedures are shown in Fig. 8for both periodic arrangements. We observe variationswith the dislocation periodic arrangement used in thesimulation, as well as with the size of the simulationcell. The results are nevertheless close to the predictionsof elasticity theory ( § III C) for all the three interactionmodels. We could even see with the EAM potential thatthe dissociation lengths extracted from atomistic simu-lations converge to the value given by elasticity theoryfor large enough unit cells. It is thus relevant to describethe screw dislocation as dissociated in two partial disloca-tions linked by a stacking fault, although the dissociationlength remains small.
C. Core energy
We obtain the dislocation core energy by subtractingthe elastic energy from the excess energy given by the D i ss o c i a ti on l e ng t h ( Å ) D i ss o c i a ti on l e ng t h ( l p = c / ) Number of atomsPwscf: SSiesta: SEAM: S OOO
FIG. 8. Dissociation length of the screw dislocation. Sym-bols correspond to data extracted from atomistic simulationsthrough the fit of the disregistry (Eq. 8) for both the O andS periodic arrangements and for different sizes of the simu-lation cell. The solid lines are the predictions of elasticitytheory based on the stacking faults (Tab. II). On the rightvertical axis, the experimental c lattice parameter has beenused to normalized the dissociation length by the distance λ P = c/
120 140 160 18010 E c o r e ( m e V Å − ) Number of atomsPwscf: SSiesta: SEAM: S OOO
FIG. 9. Core energy of the screw dislocation obtained fordifferent sizes of the atomistic simulation cell and for boththe O and S periodic arrangements ( r c = b ). atomistic simulations. This excess energy is the energydifference between the simulation cell containing the dis-location dipole and the same cell without any defect.The elastic energy calculation takes into account theelastic anisotropy and the effect of periodic bound-ary conditions. The obtained core energy are shown inFig. 9. Like for the dissociation length, some variationswith the size of the simulation cell and the periodic ar-rangement can be observed. We did not manage to linkboth quantities, i.e. the core energy and the dissociationlength. They are not simply related by the expression ofthe energy variation with the dissociation length (Eq. 3),and the change of the elastic interaction between disloca-tions caused by their dissociation could not fully explainthe variation of the core energy either. As the spread- ing of partial dislocations also depends on the size of thesimulation cells, the variation of the core energy proba-bly have a more complex origin than simply a variationof the dissociation length.We could use, with the EAM potential, a simulationcell large enough to obtain a converged value of the coreenergy, 163 meV ˚A − for a core radius r c = b . It is worthpointing that this value does not depend on the periodicarrangement used in the simulation (Fig. 9). This can beachieved thanks to a proper account of the core tractioncontribution to the elastic energy. A difference of ∼
10 meV ˚A − would have been observed between the S andO periodic arrangement without this contribution. Abinitio calculations lead to a dislocation core energy of145 ± − for Siesta and 125 ± − for Pwscf ( r c = b in both cases). D. Dissociation in the basal plane
Basal slip is observed experimentally only at high tem-perature (above 850 K) and for a higher resolved shearstress than the one needed to activate prismatic slip. Atlower temperatures, no basal slip could be observed, even when the monocrystal was oriented so as to favorbasal slip. Our atomistic simulations lead to a screwdislocation configuration dissociated in the prism plane,which clearly cannot glide easily in the basal plane. It isworth looking if another configuration, which could glidein this basal plane, also exists. [101-0][0001][12-10]
FIG. 10. Differential displacement map of the metastableconfiguration of a 1 / n = 6 (288 atoms). Crosses × correspond to the positions of the two partial dislocations, and+ to their middle, i.e. the position of the total dislocation. Using the same EAM potential, Khater and Bacon showed that a screw dislocation can also dissociate inthe basal plane. The basal configuration is obtained byintroducing in the same basal plane two partial disloca-tions of Burgers vector 1 / / γ -surface (Fig. 1). After relaxation of the atomic positions,the dislocation remains dissociated in the basal plane, ascan be seen from the corresponding differential displace-ment map (Fig. 10). A fit of the screw component of thedisregistry created by the dislocation in the basal planeleads to a dissociation length d = 6 . − than the configuration dissociated in theprism plane.We check if ab initio calculations also lead to such ametastable configuration. Starting from a screw disloca-tion dissociated in two partial dislocations in the basalplane, both Siesta and
Pwscf lead after relaxation ofthe atomic positions to the stable configuration dissoci-ated in the prism plane. This is true both with the S( n = 5) and the O ( n = 4) periodic arrangements. These ab initio calculations show then that such a dissociationof the screw dislocation in the basal plane is unstable.The metastable configuration observed with the EAMpotential appears to be an artifact of the empirical po-tential. This is not specific to the Mendelev and Acklandpotential as a configuration dissociated in the basalplane is stabilized by any central forces potential. V. PEIERLS BARRIER
Before calculating ab initio the Peierls barrier of thescrew dislocation, we use the EAM potential to assessthe validity of the method and to check the convergenceof the Peierls barrier with the size of the simulation cell.
A. Methodology
We determine the Peierls barrier for the screw dislo-cation gliding in the prism plane. This is done using aconstrained minimization between two adjacent equilib-rium configurations of the dislocations. We move bothdislocations in the same direction by one Peierls distance λ P = c/ In-termediate configurations are built by linearly interpo-lating the atomic coordinates between the initial and fi-nal states. We define the corresponding reaction coor-dinate ζ = ( ~X − ~X I ) · ( ~X F − ~X I ) / k ~X F − ~X I k , where ~X , ~X I , and ~X F are the 3 N vectors defining atomic posi-tions for respectively the intermediate, initial, and finalconfigurations. In the drag method, the minimization isperformed on all atomic coordinates with the constraintthat ζ remains fixed for each of the nine intermediateimages. Nine intermediate images are also used in NEBmethod with a spring constant k = 0 . − .So as to obtain the variation along the path of the dis-location energy with its position, we need to determinethe dislocation position x Dislo ( ζ ) for each intermediateimage, once it has been relaxed. This is done thanks toa fit of the disregistry in the prism plane, as describedin § IV B. This allows us to check that both dislocations D E P ( m e V Å − ) x Dislo / l P NEBDrag 0 0.25 0.5 0.75 0 0.25 0.5 0.75 1 x D i s l o / l P z NEBDragx
Dislo = z l P D E P ( m e V Å − ) NEBDrag
FIG. 11. Peierls barrier of a screw dislocation calculated withthe EAM potential for the S periodic arrangement with 1600atoms ( n = m = 20). Two different methods for findingthe minimum energy path have been used: the NEB methodand a simple constrained minimization (drag). Symbols cor-respond to the results of atomistic simulations and lines totheir interpolation with Fourier series. in the simulation cell are moving in a coordinated way:the distance between them remains fixed. As a conse-quence, there is no variation of the elastic energy alongthe path and the energy variation ∆ E P ( ζ ) obtained bythe constrained minimization corresponds to a variationof the core energy, i.e. the Peierls energy. We deduce thePeierls energy ∆ E P ( x Dislo ) by eliminating the reactioncoordinate ζ between ∆ E P ( ζ ) and x Dislo ( ζ ).Fig. 11 illustrates the whole procedure for a given dislo-cation periodic arrangement using both constrained min-imization techniques. The variation ∆ E P ( ζ ) slightly dif-fers between both techniques: for a given reaction coor-dinate ζ , the drag method leads to a state of lower energythan NEB method. This corresponds to a small differ-ence in the dislocation position: this position deviatesmore with drag than with NEB method from a linearvariation. Nevertheless, one obtains at the end the samePeierls barrier ∆ E P ( x Dislo ), whatever the method used.These differences observed for the functions ∆ E P ( ζ ) and x Dislo ( ζ ) between drag and NEB methods increase withthe size of the simulation cell, i.e. with the number ofdegrees of liberty. For large simulation cells (contain-ing more than 3600 atoms in the S periodic arrangementfor instance), the drag method sometimes fails to find acontinuous path between the initial and final states: onehas to use the NEB method in theses cases. Only much0 D E P ( m e V Å − ) x Dislo / l P (a) S arrangement 100 atoms14425640016003600 0 0.1 0.2 0.3 0.4 0 0.25 0.5 0.75 1 D E P ( m e V Å − ) x Dislo / l P (b) O arrangement 392 atoms8003200
FIG. 12. Variation of the Peierls barrier with the size of thesimulation cell calculated with the EAM potential for bothperiodic arrangements. smaller simulation cells can be studied ab initio . Forthese sizes, the drag and the NEB methods always leadto the same result. We will therefore only use the dragmethod in the ab initio calculations, as it costs much lessCPU time.Finally, we interpolate with Fourier series the periodicfunctions ∆ E P ( ζ ) and x Dislo ( ζ ) − λ P ζ . This leads to asmooth function ∆ E P ( x Dislo ) that can be derived. ThePeierls stress σ P is deduced from the maximal slope ofthis function, σ P = 1 b Max (cid:18) ∂ ∆ E P ∂x Dislo (cid:19) . (9)We obtain a Peierls stress σ P = 24 ± determined, for thesame empirical potential, a Peierls stress σ P = 22 MPausing molecular statics simulations under applied stress.As the agreement between both methods is good, we seethat the Peierls stress can be defined either from the slopeof the Peierls barrier or from the minimal applied stressunder which the dislocation glides indefinitely.We now examine, still with the EAM potential, howthis Peierls barrier varies when the size of the simulationcell decreases up to reaching a number of atoms that canbe handled in ab initio calculations (Fig. 12). Both the S and O periodic arrangements give the same Peierls bar-rier, and hence the same Peierls stress, for a large enoughsimulation cell ( & ab initio this Peierls barrier.This will allow us to confirm that the Peierls stress is aslow as indicated by experiments and this EAM poten-tial. Moreover, it is worth pointing that the expectedPeierls barrier should be small: ∆ E P = 0 . − atthe saddle point according to the EAM potential. Thiscorresponds to a difference of energy 2 b ∆ E P = 2 . ab initio calcula-tions. Looking for an upper limit of this value is easier asit minimizes the problems associated with this precision.Finally, it is worth pointing that the dissociation lengthvaries during the dislocation migration, and this variationis more pronounced ( ∼ B. Ab initio barriers
The Peierls barriers obtained by ab initio calculationsare shown on Fig. 13(a) for
Pwscf and Fig. 13(b) for
Siesta . In both cases, the height of the barrier decreaseswhen the number of atoms increases, in agreement withwhat has been observed with the EAM potential for thesame S periodic arrangement. For a given number ofatoms, the ab initio barriers are a little bit lower than theones obtained with the EAM potential. Results obtainedwith
Siesta are noisy: the energy barrier that we wantto calculate is so small that it needs a really strict conver-gence criterion on atomic forces for the relaxation. Wedid not manage to reach such a precision with
Siesta .On the other hand, the barriers obtained with
Pwscf are smooth. We can therefore interpolate the ab initio results with Fourier series, and estimate then the Peierlsstress (Eq. 9). This leads to σ P = 36 MPa for the simula-tion cell containing 100 atoms and σ P = 21 MPa for 144atoms. Considering that these values, obtained in smallsimulation cells, are upper limits, ab initio calculationspredict a Peierls stress smaller than 21 MPa for a screwdislocation in zirconium gliding in a prism plane.The comparison of this ab initio estimate of the Peierlsstress with experimental data is quite challeng-ing. As pointed out in the introduction, the yield stressof zirconium strongly depends on its oxygen content.No experimental data exists for a purity better than0.07% O (in atomic fraction). Moreover, all measure-1 D E P ( m e V Å − ) x Dislo / l P (a) Pwscf 100 atoms144 0 0.2 0.4 0.6 0.8 0 0.25 0.5 0.75 1 D E P ( m e V Å − ) x Dislo / l P (b) Siesta 100 atoms144 atoms196 atoms
FIG. 13. Peierls barrier for a screw dislocation gliding inits prism plane calculated ab initio with (a)
Pwscf and (b)
Siesta for the S periodic arrangement. Symbols correspondto ab initio results and lines to their interpolation by Fourierseries. t ( M P a ) O content (at. %)Rapperport and Hartley (1960)Mills and Craig (1968)Soo and Higghing (1968)Akhtar and Teghtsoonian (1971)
Ab initio upper limit
FIG. 14. Zirconium flow stress determined experimentally forvarious O contents.
Data have been extrapolated to0 K and are compared to the ab initio
Peierls stress of screwdislocations in pure Zr. ments have been performed at temperatures higher than77 K. The determination of pure zirconium flow stress at0 K therefore needs some extrapolation of experimentaldata. Without a clear understanding of the mechanismsresponsible of zirconium hardening by oxygen impurities,such an extrapolation is quite hazardous. Nevertheless, agraphical comparison (Fig. 14) of experimental data withour ab initio
Peierls stress shows a reasonable agreement,keeping in mind that the value 21 MPa has to be consid-ered as an upper limit.
VI. CONCLUSIONS
Using two ab initio approaches (
Pwscf and
Siesta ),both in the DFT-GGA approximation, we have shownthat a 1 / h i screw dislocation in hcp zirconium dis-sociates in two partial dislocations with a pure screwcharacter. This is in agreement with the minimum in1 / h i found for the generalized stacking fault energyin the prism plane. We could extract the dissociationlength from our atomistic simulations. Although this dis-sociation length is small ( d ∼ Pwscf and d ∼ Siesta ), it is in reasonable agreement with the onepredicted by elasticity theory.The EAM potential of Mendelev and Ackland leadsto the same structure of the dislocation dissociated in theprism plane. The metastable configuration dissociated inthe basal plane predicted by this potential is not stablein ab initio calculations, both with Pwscf and
Siesta .This configuration is therefore an artifact of this poten-tial, probably induced by the deep minimum in 1 / h i found with this empirical potential for the generalizedstacking fault in the basal plane. This minimum is muchmore shallow in ab initio calculations.We could also obtain an ab initio estimate of thePeierls stress of the screw dislocation gliding in the prismplane. Calculations with Pwscf lead to an upper limitof 21 MPa for this Peierls stress. This small value showsthat screw dislocations can glide quite easily in pure zir-conium, thus confirming what had been obtained withMendelev and Ackland EAM potential. Such a smallPeierls stress is in agreement with experimental data,once the hardening of oxygen impurities has been con-sidered.
ACKNOWLEDGMENTS
This work was performed using HPC resources fromGENCI-CINES and GENCI-CCRT (Grant Nos 2011-096020 and 2012-096847). ∗ [email protected] E. Rapperport, Acta Metall. , 254 (1959). P. Soo and G. T. Higgins, Acta Metall. , 177 (1968). A. Akhtar and A. Teghtsoonian, Acta Metall. , 655 (1971). D. Caillard and J. L. Martin,
Thermally Activated Mechanisms in Crystal Plasticity ,edited by R. W. Cahn (Pergamon, Amsterdam, 2003). B. Legrand, Philos. Mag. B , 171 (1984). B. Legrand, Philos. Mag. A , 43 (1986). B. Legrand, Philos. Mag. A , 83 (1985). A. Girshick, D. G. Pettifor, and V. Vitek,Philos. Mag. A , 999 (1998). F. Ferrer, A. Barbu, T. Bretheau, J. Cre-pin, F. Willaime, and D. Charquet, in
Zirconium in the nuclear industry: thirteenth international symposium ,American Society for Testing and Materials Special Tech-nical Publication, Vol. 1423, edited by G. D. Moan andP. Rudling (American Society Testing and Materials, WConshohocken, USA, 2002) pp. 863–885, 13th Interna-tional Symposium on Zirconium in the Nuclear Industry,Annecy, France, June 10-14, 2001. C. Domain and A. Legris, in
Mat. Res. Soc. Symp. Proc. ,Vol. 653 (2001) p. Z3.8. C. Domain, J. Nucl. Mater. , 1 (2006). N. Tarrat, M. Benoit, and J. Morillo,Int. J. Mater. Res. , 329 (2009). M. Ghazisaeidi and D. Trinkle,Acta Mater. , 1287 (2012). V. Vitek and V. Paidar, in
Dislocations in Solids , Vol. 14,edited by J. Hirth (Elsevier, 2008) Chap. 87, pp. 439–514. J. Bailey, J. Nucl. Mater. , 300 (1962). S. Naka, A. Lasalmonie, P. Costa, and L. P. Kubin,Philos. Mag. A , 717 (1988). S. Farenc, D. Caillard, and A. Couret,Acta Metall. Mater. , 2701 (1993). S. Farenc, D. Caillard, and A. Couret,Acta Metall. Mater. , 3669 (1995). D. Mills and G. B. Craig, Trans. AIME , 1881 (1968). D. H. Sastry, Y. V. R. K. Prasad, and K. I. Vasu,J. Mater. Sci. , 332 (1971). D. Sastry and K. Vasu, Acta Metall. , 399 (1972). T. Tanaka and H. Conrad, Acta Metall. , 1019 (1972). A. Akhtar and E. Teghtsoonian,Metall. Mater. Trans. A , 2201 (1975). M. I. Mendelev and G. J. Ackland,Philos. Mag. Lett. , 349 (2007). H. A. Khater and D. J. Bacon,Acta Mater. , 2978 (2010). J. M. Soler, E. Artacho, J. D. Gale, A. Garc´ıa,J. Junquera, P. Ordej´on, and D. S´anchez-Portal,J. Phys.: Condens. Matter , 2745 (2002). P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fab-ris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougous-sis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcov-itch, J. Phys.: Condens. Matter , 395502 (2009). C. Domain, R. Besson, and A. Legris,Acta Mater. , 1495 (2004). G. V´erit´e, F. Willaime, and C. C. Fu,Solid State Phenom. , 75 (2007). C. S. G. Cousins, J. Phys. C Solid State , 989 (1979). P. Villars and L. D. Calvert,
Pearson’s Handbook of Crys-tallographic Data for Intermetallic Phases (American So-ciety for Metals, OH, 1985). C. Stassis, J. Zarestky, D. Arch, O. D. McMasters, andB. N. Harmon, Phys. Rev. B , 2632 (1978). E. S. Fisher and C. J. Renken,Phys. Rev. , A482 (1964). Elastic constants calculated by Mendelev and Ackland inthe original article describing the EAM potential did nottake into account atomic relaxations. This explains the dif-ference with the ones given in Tab. I. V. Vitek, Philos. Mag. , 773 (1968). Y. Udagawa, M. Yamaguchi, H. Abe, N. Sekimura, andT. Fuketa, Acta Mater. , 3927 (2010). D. Hull and D. J. Bacon,
Introduction to Dislocations , 4thed. (Butterworth Heinemann, 2001). See Supplemental Material athttp://link.aps.org/supplemental/10.1103/PhysRevB.86.144104for a discussion on the convergence of this fault energywith
Pwscf parameters. G. J. Ackland, S. J. Wooding, and D. J. Bacon,Philos. Mag. A , 553 (1995). D. Bacon and V. Vitek,Metall. Mater. Trans. A , 721 (2002). K. Schwartzkopff, Acta Metall. , 345 (1969). J. P. Hirth and J. Lothe,
Theory of Dislocations , 2nd ed.(Wiley, New York, 1982). A. N. Stroh, Philos. Mag. , 625 (1958). A. N. Stroh, J. Math. Phys. (Cambridge, Mass.) , 77(1962). A. J. E. Foreman, Acta Metall. , 322 (1955). L. J. Teutonico, Mater. Sci. Eng. , 27 (1970). M. M. Savin, V. M. Chernov, and A. M. Strokova,Phys. Status Solidi A , 747 (1976). Foreman gave a different expression for K , but thecomparison with a numerical evaluation of the K matrixusing Stroh formalism shows that the correct expres-sion is the one given by Savin et al .. . V. V. Bulatov and W. Cai,
Computer Simulations of Dis-locations , edited by A. P. Sutton and R. E. Rudd, Oxfordseries on materials modelling (Oxford University Press,2006). W. Cai, V. V. Bulatov, J. Chang, J. Li, and S. Yip,Philos. Mag. , 539 (2003). V. Vitek, R. C. Perrin, and D. K. Bowen,Philos. Mag. , 1049 (1970). R. Peierls, Proc. Phys. Soc. , 34 (1940). F. R. N. Nabarro, Proc. Phys. Soc. , 256 (1947). E. Clouet, Philos. Mag. , 1565 (2009). A. Akhtar, Acta Metall. , 1 (1973). K. E. J. Rapperport and C. S. Hartley, Trans. AIME ,869 (1960). M. H. Liang and D. J. Bacon,Philos. Mag. A , 181 (1986). G. Henkelman, G. J´ohannesson, and H. J´onsson,in