Calculation of dislocation positions and curved transition pathways in BCC crystals from atomic displacements
aa r X i v : . [ c ond - m a t . m t r l - s c i ] S e p Calculation of dislocation positions and curved transitionpathways in BCC crystals from atomic displacements
R. Gr¨oger ∗ Central European Institute of Technology – Institute of Physics of Materials (CEITEC-IPM),Academy of Sciences of the Czech Republic, ˇZiˇzkova 22, 61662 Brno, Czech Republic (Dated: April 7, 2019)The thermodynamic description of dislocation glide in crystals depends crucially on the shapeof the Peierls barrier that the dislocation has to overcome when moving in the lattice. While theheight of this barrier can be obtained unequivocally using saddle-point search algorithms such asthe Nudged Elastic Band (NEB) method, its exact shape depends on the chosen approximationof the transition pathway of the system. The purpose of this paper is to formulate a procedurethat allows to identify the position of the dislocation directly from the displacements of atoms inits core. We investigate the performance of this model by calculating curved paths of a 1 / PACS numbers: 61.72.Hh, 02.60.-x, 45.10.DbKeywords: screw dislocation, BCC metal, dislocation pathway, Peierls barrier, Nudged Elastic Band.
I. INTRODUCTION
Computational studies of dislocations in body-centeredcubic (BCC) metals have been made using a wide vari-ety of interatomic potentials and ab initio methods .They view the evolution of the state of the system as itsmotion along an a priori unknown path in the configura-tional space spanned by the 3 N atomic degrees of free-dom (DOF), where N is the number of atoms in the simu-lated block. While this picture arises directly from molec-ular simulations, it presents severe difficulties when de-veloping theoretical models of thermally activated glideof dislocations , which serve to coarse-grain the atom-istic results to the continuum. Not surprisingly, there hasbeen a long-standing interest in developing approximateschemes to extract the effective position of the dislocationfrom the positions of atoms obtained by molecular staticssimulations of single dislocations. While this approach iscertainly attractive, it cannot be used without developinga systematic procedure that maps the 3 N atomic DOFto the position of the dislocation.The most obvious choice that leads to the reduction ofcomplexity of the system is to assume that the dislocationmoves between two neighboring lattice sites along thestraight line. This is implicitly assumed in most papersemploying the Nudged Elastic Band (NEB) method ,where the position of the dislocation along the minimumenergy path scales linearly with the image number .While the obtained barrier can be used to assess the sta-bility of each intermediate state along the path, it doesnot constitute the Peierls barrier that could be used todevelop thermodynamic models of dislocation glide suchas that due to Dorn and Rajnak . Moreover, a straightapplication of the NEB method to all DOF in the sys-tem leads to nonuniform distributions of dislocation po-sitions among the images, which affects the shape of the obtained Peierls barrier . A significant improvement ofthese results is obtained using our modification of theNEB method , where atomic relaxations are taken intoaccount. This NEB+r method guarantees that the po-sitions of the dislocation when following the minimumenergy path are distributed uniformly among the imagesand thus the assumed proportionality between the dislo-cation position and the image number is justified. Thesedevelopments have led to an accurate estimate of thePeierls barrier of 1 / h i screw dislocations in BCC Wand its changes under the applied stress . However,these calculations are still based on the assumption thatthe path of the dislocation between two neighboring im-ages is a straight line, which is not true in general.In principle, it should be possible to deduce the dislo-cation position (and thus also its path between two po-sitions in the lattice) directly from the displacements ofatoms as obtained from molecular statics calculations orfrom the NEB (NEB+r) methods. This idea dates backto Peierls and Nabarro , who associated the disloca-tion position with the point in the slip plane at whichthe displacement parallel to the slip direction, interpo-lated from the displacements of atoms, is equal to b/ b is the magnitude of the Burgers vector of thedislocation. Since then, this argument was used manytimes. In particular, it was adopted by Pizzagalli et al. ,Rodney and Proville and Proville et al. , where theposition of the dislocation is defined by a single coordi-nate corresponding to the distance that the dislocationmakes in a well-defined slip plane. The same level of ap-proximation was also used by Ventelon et al. in one oftheir methods (disregistry function method) that uses acombination of isotropic elasticity and geometry of theBCC lattice to define the position of the dislocation.In general, the movement of the dislocation should beviewed as a three-dimensional event during which thecenter of the dislocation transits along a curved pathin the vicinity of the slip plane. In attempt to resolvethis path, Ventelon et al. also proposed another model,whereby the dislocation position is identified using a costfunction that is based on both the actual positions ofatoms (as obtained from atomistic simulations) and thepositions of atoms obtained from anisotropic elasticityfor some trial dislocation position in the { } plane.The actual position of the dislocation is then obtainedby minimizing this cost function that is defined as thedistance between the two sets of coordinates in the five-dimensional subspace spanned by the coordinates of thefive most displaced atoms around the dislocation. Sincethe position of the dislocation is determined by relativedisplacements of the three atoms closest to the disloca-tion in the direction parallel to the Burgers vector, a sim-ilar approach can be devised that is based on inversionof the Eshelby-Stroh sextic formalism (see, for exampleHirth and Lothe ), which provides elastic displacementsof atoms corresponding to a given position of the dis-location. We have investigated this possibility earlier(unpublished work). While this approximation can beused when the dislocation is in the middle of the latticesite formed by the nearest three h i atomic columns,it quickly worsens as the dislocation gets closer to any ofthese columns or the boundary between the neighboringlattice sites.A different scheme whereby the position of the disloca-tion is determined by extrapolating differential displace-ments between the three atoms surrounding the dislo-cation into the interior of this triangle was developedby Itakura et al. . This approach is closely relatedto a purely geometrical concept of barycentric (or tri-linear) coordinates known from ternary diagrams thatwas originally applied to estimate the position of the dis-location by Heinrich and Schellenberger . While thismethod makes use of the actual positions of atoms, theexpression of the dislocation position as a linear combi-nation of the displacements of the three nearest atoms inthe direction parallel to the Burgers vector represents aconvenient choice that is, however, not justified phys-ically. This approximation was avoided in the recentwork of Dezerald et al. who aimed to reconstruct thetwo-dimensional Peierls barrier by interpolating the lineenergy of the dislocation, calculated by first principles,from two straight dislocation paths. One connects theneighboring potential minima in the { } plane and theother passes from an atom (“split-core” configuration) tothe so-called “hard-core” position in the h i directionperpendicular to the first path. The Peierls potential isthen expressed in the form of a Fourier series with theFourier coefficients adjusted so as to minimize the leastsquares error between the calculated data and the Fourierseries. The path of the dislocation between two neighbor-ing minimum-energy lattice sites in the { } plane areobtained using the disregistry and cost function methodsof Ventelon et al. . While these paths are smooth, thecorresponding Peierls barriers display sharp maxima for BCC Mo, W, Nb, while this is somewhat less pronouncedfor BCC Ta, V and ferromagnetic BCC Fe. This does notagree with the work of Suzuki et al. and our more recentwork (Ref. 13,27), which show that in order to reproducethe experimentally measured temperature dependence ofthe flow stress, the Peierls potential has to possess a flatmaximum.In this paper, we develop a procedure that providesthe position of a 1 / in that it considers all three { } planes onwhich the dislocation can move. These calculations pro-vide three lines whose intersection defines the positionof the dislocation in the perpendicular (111) plane. Wedemonstrate the performance of this method by calcu-lating the paths of a straight 1 / { } planes from thediscrete snapshots (images) of the system obtained re-cently using NEB+r calculations . These calculationsare made under zero applied stress and for positive andnegative shear stresses perpendicular to the slip direc-tion. We demonstrate that the shape of the Peierls bar-rier changes when considering the curved transition path-way of the dislocation as compared to that obtained pre-viously by assuming the straight path of the dislocationbetween two equivalent minimum-energy configurationsin the lattice . II. COMPLEXITY REDUCTION
Let us consider that X i is the position of atom i in the perfect lattice. Upon inserting a 1 / Σ and relaxing the atoms,each atom i moves into its new position, x i , where i = 1 , , . . . , N are atomic numbers. This configurationcorresponds to the stable equilibrium (a minimum) withenergy E ( x . . . x N ; Σ ). In this configuration the dislo-cation is at the bottom of the Peierls valley for the ap-plied stress Σ . We will now consider that the dislocationhas moved away from this minimum at constant appliedstress Σ , which is to say that all atoms were displacedfrom x i to x i . The energy of this new (nonequilibrium)configuration will be denoted E ( x . . . x N ; Σ ). The en-thalpy of the final state of the system relative to its initialstate per unit length of the dislocation (or, simply, theline enthalpy of the dislocation) is then H ( x . . . x N ; Σ ) = E ( x . . . x N ; Σ ) − E ( x . . . x N ; Σ ) l dislo , (1)where l dislo is the length of the dislocation segment con-tained in the simulated block.The left-hand side of (1) seems to imply that the stateof the system is described by the 3 N DOF associated withthe positions of all particles. However, this is rather im-practical from the computational point of view because,in this case, the energy of the system would have to becalculated from all atoms in the system.To reduce this complexity, it is more convenient to sep-arate the atomic degrees of freedom into two classes: (i)the minimum number of DOF that determine the posi-tion of the dislocation, and (ii) all remaining DOF thatdo not affect the position of the dislocation and thus theyrepresent merely the large-scale response of the systemto incorporating the dislocation. For this purpose, wecan use the concept of the Burgers circuit to identifythe region of the atomic block that contains the disloca-tion in its interior. When viewing a BCC crystal alongthe [111] direction, the shortest of these circuits passesthrough three atoms that are closest to the center of thedislocation. This implies that only 9 DOF are necessaryto describe the dislocation position, while the remaining3 N − X D , Y D ) isencoded in only two DOF.The former suggests to find a procedure that mapsthe positions of atoms, x . . . x N , to the position of theintersection of the dislocation line with the perpendicular(111) plane in the perfect lattice (called hereafter as thedislocation position), i.e. M : { x . . . x N } 7→ ( X D , Y D ) . (2)It is important to emphasize that the former is defined inthe deformed lattice with the dislocation, while the latterin the perfect lattice. If this mapping exists, the state ofthe system can be described using only two variables,the coordinates ( X D , Y D ) of the dislocation in the (111)plane of the perfect lattice. Hence, one can write theline enthalpy of the dislocation as H ( X D , Y D ; Σ ) = H ( x . . . x N ; Σ ) , (3)where the right-hand side is obtained from (1). Eq. (3)opens the possibility for developing a model of thermallyactivated dislocation glide, in which the transformationof the dislocation core is viewed as a motion of the centerof the dislocation in the underlying Peierls potential.For further developments, it will be convenient to in-troduce a transition pathway of the dislocation, ξ , that isdefined by a series of discrete points ( X D , Y D ) obtainedfrom individual snapshots of the system as it moves alongthe minimum energy path. Hence, the left-hand side of(3) is equivalent to H ( ξ ; Σ ), H ( ξ ; Σ ) ≡ H ( X D , Y D ; Σ ) , (4) where ξ represents a particular point ( X D , Y D ) along acurved transition path of the dislocation. If the disloca-tion remains a straight line during this transition (whichis the case at 0 K), we may express the line enthalpy ofthe dislocation as H ( ξ ; Σ ) = V ( ξ ; Σ nonglide ) − σ glide bξ , (5)where V is the Peierls barrier and the second term is thework done by the applied stress on displacing the disloca-tion by the distance ξ measured along the transition path.In the equation above, we have split the applied stress as Σ = Σ glide + Σ nonglide , where Σ glide contains only theshear stress σ glide acting in the slip plane parallel to theslip direction (i.e. the Schmid stress), whereas Σ nonglide contains all other stresses. The latter are all stress com-ponents that do not exert a Peach-Koehler force on thedislocation. Clearly, the glide (Schmid) stress does workon displacing the dislocation, while non-glide stresses af-fect the shape of the Peierls barrier. The Peierls barriercan then be obtained from V ( ξ ; Σ nonglide ) = H ( ξ ; Σ ) + σ glide bξ . (6)Here, the first term on the right-hand side is obtaineddirectly from the NEB (NEB+r) calculations, as it isevident by combining (1), (3) and (4). The expression(6) is completely general and can be used to obtain theshape of the Peierls barrier for an arbitrary applied load.For the sake of simplicity, we will consider in the fol-lowing that σ glide = 0. Hence, the NEB (NEB+r) cal-culations directly yield the line energies of the disloca-tion, V ( ξ ; Σ nonglide ), subject to a given non-glide stress Σ nonglide . III. POSITION OF THE DISLOCATION
Upon inserting the dislocation parallel to the z axis,applying external load and relaxing the atomic positions,the atoms move from their perfect lattice positions, X i ,into their new positions, x i . These two sets of atomic po-sitions can be used to obtain a differential displacementmap that uniquely identifies the lattice site with the dis-location. One such site corresponds to the gray trianglein Fig. 1 with the dislocation at some unknown positionin this interior. The three { } planes on which the dis-location can move are marked in Fig. 1 as α = 1 , , α , we consider the two planes p ± α immediately above and below the lattice site with thedislocation, which are drawn in Fig. 1 by solid lines. Each of these planes is represented by a finite number n ± α of atoms P ± α ( i ) , where i = 1 . . . n ± α , which are orderedsuch as to form chains. The displacement of each atomprojected parallel to the Burgers vector of the dislocationis then defined as u b ( P ± α ( i ) ) = h x ( P ± α ( i ) ) − X ( P ± α ( i ) ) i · b b , (7) FIG. 1: The three possible { } slip planes on which thedislocation can move in the lattice ( α = 1 , , p ± α . The dislocationpositions within these chains are marked Q ± α . The sought po-sition of the dislocation corresponds to the point D , which isdefined as an intersection of the three line segments Q + α Q − α for α = 1 , , where b is the Burgers vector of the dislocation, and b itsmagnitude. For the atoms far away from the dislocationcore, the displacements u b approach constant values thatrepresent the maximum and minimum of u b along thechain. Similarly as in Refs. 18,19, the position of thedislocation in each chain p ± α , denoted hereafter Q ± α , isassumed to coincide with the point along the chain, wherethe displacement in the direction parallel to the Burgersvector is half-way between its minima and maxima: u b ( Q ± α ) = 12 h min i u b ( P ± α ( i ) ) + max i u b ( P ± α ( i ) ) i . (8)For each slip plane α , this equation thus provides twoestimates of the dislocation position, one for the chainof atoms above ( Q + α ) and the other below ( Q − α ) the lat-tice site with the dislocation. These points then definethe line segments Q + α Q − α that are plotted in Fig. 1 bydashed lines. The actual position of the dislocation isassociated with the intersection of these line segments.This is marked in Fig. 1 as the point D . The implicitdefinition of Q ± α using Eq. (8) is an approximation that,nevertheless, provides a good estimate of the position ofthe dislocation when it is at the body center of the shadeddislocation triangle in Fig. 1, i.e. when the system is inequilibrium at zero applied stress. This way of obtain-ing Q ± α is adopted here for all states of the system alongthe minimum energy path irrespective of whether theyrepresent equilibrium or nonequilibrium states.It should be pointed out that there is no a priori guar-antee that the three line segments Q + α Q − α for α = 1 , , (a)(b)FIG. 2: Approximation of the dislocation position from theimage obtained using the NEB+r method, where the disloca-tion is very close to the boundary between two neighboringlattice sites (a). The small triangles are estimates of the dis-location position obtained using the method described here.The vertex marked “false node” [red circle in (a)] lies out-side of the gray triangle and represents a false prediction ofthe dislocation position, as described in the text. In (b), weshow the approximation of the dislocation position along thetransition path. The black curve marked “true nodes” rep-resents the correct approximation of the dislocation pathway.For comparison, we plot in green (+ symbols) the path inBCC W obtained by Dezerald et al. using their cost func-tion method; this is taken from their Fig.10(b) and plotted inthe orientation of the block shown in Fig. 1. three points that correspond to the intersections of threepairs of lines Q + α Q − α × Q + β Q − β for the slip planes α = β .This is illustrated in Fig. 2(a), where the shaded spherescorrespond to the positions of atoms in the perfect lat-tice. The arrows are relative displacements of atoms be-tween the image along the path for which the dislocationis closest to the boundary between the two lattice sites inFig. 2(a) and the perfect lattice. The three edges of thesmall triangles are obtained from the method describedabove. In particular, the black edge in Fig. 2(a) lies along Q +1 Q − , purple along Q +2 Q − and blue along Q +3 Q − linesdefined in Fig. 1. Each vertex of this small triangle, aswell as any point in its interior, that are within the largegray triangle can be taken as representatives of the dis-location position. This is, however, not the case for thevertex marked “false node”, because it is outside of thegray lattice site determined by the differential displace-ment map. Fig. 2(b) shows the transition pathway pre-dicted by the “true nodes” (black dots) and “false nodes”(red dots) from a series of images obtained by the NEB+runder zero applied stress . The two sets of nodes coin-cide when the system is close to the beginning, end andthe middle of the minimum energy path. For these im-ages, the three line segments in Fig. 1 intersect at a singlepoint and thus the dislocation positions are determineduniquely.We define the transition pathway as a curve connect-ing the “true nodes” in Fig. 2(b). This is quite differ-ent from that obtained by the cost function method em-ployed in Ref. 25 that is plotted by the green line (+symbols). This deviation is the largest for the point inthe middle of the path, which is determined uniquely byour method. Dezerald et al. also show that the paths ofthe dislocation at zero applied stress obtained for variousBCC metals are quite similar. This opens the possibilityfor a qualitative comparison of the dislocation pathwayshown in Fig. 2(b) with Fig. 11 (data marked “DFT”)and Fig. 14 of Itakura et al. for BCC Fe. The lat-ter show that the energy of the dislocation is the lowestwhen positioned between the points marked H (“hard-core” position) and M in Fig. 2(b), which agrees wellwith the path predicted by our method.The origin of the “false node” in Fig. 2(a) can be un-derstood by looking at the displacements u b ( P ± ) of thechain of atoms that are parallel to the boundary be-tween these neighboring lattice sites (i.e. the blue linesin Fig. 1). These displacements are shown in Fig. 3for two principally different atomic configurations. Thefirst, marked “equilibrium state”, corresponds to the ini-tial image along the path, where the dislocation is ex-actly at the body center of the gray triangle in Fig. 2(a).The second, marked “nonequilibrium state”, correspondsto the image in which the dislocation is closest to theboundary with the neighboring lattice site [case shownin Fig. 2(a)]. In the former, the displacements u b ( P ± )vary monotonously from their minimum at one end ofthe chain to their maximum at the other. In this case,the average displacement of atoms parallel to the Burg-ers vector is close to the inflection point of u b ( P ± ). Pro-vided the system is close to equilibrium, the position ofthe dislocation along the chain can thus be safely deter-mined using (8). For the nonequilibrium state, as thedislocation gets closer to the boundary between the twoneighboring lattice sites, one atomic bond in the chainthat coincides with this boundary becomes stretched toaround b/
2. This is shown in Fig. 3 and this relative dis-placement of atoms is the bigger the closer the disloca-tion is to the boundary. In this case, the displacementsof atoms u b ( P − ) do not vary monotonously along thechain. Instead, the sudden increase of displacements inthe vicinity of the dislocation is followed by a gradual de-cay of the displacements, eventually reaching a constantvalue.Because the position of the inflection point of u b ( P − )cannot be determined with sufficient accuracy fornonequilibrium states, the linear approximation (8) maybe a poor representation of the position Q − . Due to thenonmonotonous character of u b ( P − ) it is also not ad-visable to look for higher order schemes. Instead, theseobservations suggest to make the prediction of the dislo-cation position only from two { } planes, excluding theone that is parallel to the boundary between the adjacentlattice sites (in Fig. 2, this is the case for α = 3 corre-sponding to the blue edge of the small triangle). We have FIG. 3: Displacements of atoms parallel to the Burgers vec-tor of the dislocation in the chains p − and p +3 shown in Fig. 1(blue). The dashed curves (“equilibrium state”) correspondto the atomic configuration in which the dislocation is at thebody center of the lattice site (initial image for the NEBcalculation at zero applied stress), whereas the solid curves(“nonequilibrium state”) are for the image in which the dislo-cation is close to the boundary between the two neighboringlattice sites [see also Fig. 2(a)]. also checked that the jump in u b vs. P − does not appearif the dislocation position is obtained from the chain ofatoms farther away from the lattice site with the dislo-cation. However, this leads to less accurate estimates ofthe dislocation position (larger triangle obtained by theintersection of Q + α Q − α with Q + β Q − β for α = β ) and alsoto the drift of the dislocation position from the predic-tion made using the nearest chains of atoms in the three { } planes.In all calculations made in this paper, the interactionsbetween atoms were described using the Bond Order Po-tential (BOP) for BCC tungsten . However, we haveconfirmed that the nonmonotonous character of the u b vs. P − curve in Fig. 3 is not specific to the potential used.This was found by calculating the dislocation pathwayat zero applied stress with atomic interactions describedby the Ackland-Thetford potential for BCC Ta. Thisjump is also not caused by the mismatch betweeen thepositions of the far-field atoms that are held fixed duringthe NEB (NEB+r) calculations and the actual positionof the dislocation . This was verified by doubling theradius of the cylindrical block used for the NEB simu-lation, which did not produce any noticeable change inthe dependence of u b on P ± . The jump is present inthe data obtained from both NEB and NEB+r methodsand in all atomic configurations, where the dislocation isclose to the boundary between two adjacent lattice sites.It is also important to emphasize that the NEB meth-ods do not impose any constraint to keep the Burgersvector of the dislocation constant during the search forthe minimum energy path of the system. Owing to theabove-mentioned misfit of atomic positions around thedislocation and in the far-field and rounding errors, itis theoretically possible that the relative displacementsbetween the three atoms closest to the dislocation as ob-tained from the intermediate images obtained by NEB(NEB+r) calculations do not sum exactly to b . We havechecked this possibility for a few configurations, wherethe dislocation is very close to the edge of its lattice site.However, no significant deviation from the Burgers vectorbeing exactly b was detected.Finally, an important insight is gained by calculat-ing the dislocation positions when applying the stressdirectly in molecular statics calculations. In this case,pure shear stress parallel to the slip direction acting inthe (¯101) plane was applied in steps up to the value justbefore the dislocation moved (for details, see Ref. 9).For each relaxed atomic configuration at different ap-plied stresses, we examined the dependence of u b vs. P − .Interestingly, all curves were smooth without any jumpas observed from the NEB (NEB+r) calculations. Ob-viously, all atomic configurations obtained from directmolecular statics calculations are equilibrium states ofthe system at the given applied stress, while some inter-mediate images obtained using the NEB methods corre-spond to the unstable branch of the Peierls potential onwhich the dislocation moves spontaneously to its nearestminimum-energy configuration. As far as we can judgefrom our NEB+r calculations using 32 movable images,the origin of the jump in the u b vs. P − curve is asso-ciated with the inflection point of the Peierls barrier atwhich the dislocation can no longer be stabilized by theapplied stress. IV. NUMERICAL SIMULATIONS
In our previous publication (Ref. 17), we used theNEB+r method to investigate how the Peierls barrierof a 1 / { } planes, and (iii) plot the energies of thesesnapshots along these curved paths to obtain the Peierlsbarriers.We will first investigate how the shape of the disloca-tion pathway is affected by the shear stress perpendicu-lar to the slip direction, τ . We carry out these calcula-tions for three representative values of this stress, namely τ /C = {− . , , . } that is applied by imposing auniform stress tensor Σ ≡ Σ nonglide = diag( − τ, τ,
0) inthe coordinate system shown in Fig. 1 (for details, seeRef. 17). The calculated paths of the dislocation mov-ing on the three { } planes in the [111] zone and thethree values of the shear stresses perpendicular to theslip direction are shown in Fig. 4. The large trianglesrepresent the boundaries of four lattice sites – one fromwhich the dislocation makes the jump (shaded) and the three others representing the target sites into which thedislocation moves by the glide on the three { } planes.The dislocation paths in bold are obtained as piecewiselinear interpolations of the dislocation positions obtainedfor each image by the method developed above. (a) (b) (c)FIG. 4: Dislocation pathways calculated from the positionsof atoms in 15 movable images obtained using the NEB+rmethod. The three figures correspond to different shearstresses perpendicular to the slip direction: (a) τ /C = − .
04, (b) τ /C = 0, (c) τ /C = 0 .
04. The calculateddislocation paths are shown by solid lines. The colors distin-guish the paths of the dislocation on the three { } planes:(¯101) = black, (0¯11) = blue, and (¯110) = purple. The arrowsin the middle panel show the curvilinear coordinates ξ for thethree paths. The Peierls barriers for the three { } planes andthe three values of the shear stress perpendicular tothe slip direction ( τ ) are shown in the upper panels ofFig. 5. Here, dashed curves correspond to the barriersobtained in Ref. 17, where we considered that the dislo-cation moves between the two minima along the straightpath. In this case, the transition coordinate is definedas ξ = a I/ ( M + 1), where I = { , , . . . , M + 1 } areimage numbers and a = a p / { } planes( a = 3 . V ( ξ ; Σ nonglide ), obtained from (6), along the pathsidentified using the procedure developed in Section IIIand shown in Fig. 4. Firstly, the three barriers for τ = 0are identical, as dictated by symmetry, which is an im-portant test of our procedure to calculate the dislocationpathway. The shape of the Peierls barrier changes signif-icantly when drawn along the curved path. In particular,the panels for zero and positive τ show that the barriersin black that correspond to the glide of the dislocation onthe (¯101) plane are wider, with steeper gradient close to ξ = 0. The Peierls barrier thus deviates from the sin ξ form to a parabolic or flat-top potentials that are oftenused to construct models of thermally activated disloca-tion glide (for a review, see the book of Caillard andMartin ). The right panel in Fig. 5 shows that the (0¯11)and (¯110) barriers develop lower minima than the initialand final configurations along the path. These are prob-ably consequences of the fact that the initial and finalatomic positions are relaxed down to the maximum forceon atom of only 0 .
005 eV/˚A. L i n e e n e r gy V [ m e V / Å ] τ / C =-0.04 τ / C =0 τ / C =+0.040 1 2 ξ [Å]-100-50050100 d V / d ξ [ m e V / Å ] ξ [Å] 0 1 2 3 ξ [Å]CRSS b CRSS b /2(101) slip(110) slip(011) slip FIG. 5: Comparison of the Peierls barriers of the 1 / . The barriers plotted by dashed lines are obtained by assuming that the dislocation movesalong the straight line connecting the two neighboring potential minima in the glide plane. The solid curves are obtained byconsidering the curved path of the dislocation, as determined by the method developed in this paper and shown in Fig. 4.For τ = 0 (middle panel), the Peierls barriers corresponding to dislocation glide on the three { } planes are identical. Thebottom panel shows the derivatives of the Peierls barriers for the respective slip planes. For completeness, the lower panels of Fig. 5 show com-parisons of the derivatives of these Peierls barriers. Here,we consider only those barriers that satisfy the funda-mental equation CRSS b = max(d V / d ξ ), i.e. the dislo-cation moves on the (¯101) plane (black curves, zero orpositive τ ), or CRSS b/ V / d ξ ) if it moves onthe (¯110) plane (purple curves, negative τ ). The curva-ture of the path shifts the peak of max(d V / d ξ ) towardsthe beginning of the path. Qualitatively speaking, thisshift together with the steep increase of the Peierls bar-rier close to ξ = 0 will affect the activation enthalpy forthe dislocation glide . This can be demonstrated eas-ily in the limit of zero applied stress for which the ac-tivated segment of the dislocation contains two isolatedkinks between ξ = 0 and ξ max (the length of the curvedtransition path of the dislocation). The energy of thisactivated state relative to that of the straight dislocationat the bottom of the Peierls valley is obtained from theDorn-Rajnak model [Eq. (6) in Ref. 13 with σ ∗ = 0 and V ( ξ ) = V (0) = 0] and reads H b ≡ H k = 2 Z ξ max p [ V ( ξ )] + 2 EV ( ξ ) d ξ , (9)where E = µb / ≈ .
147 eV/˚A is the line tension of astraight dislocation . We now use the two Peierls bar-riers shown in the middle panel of Fig. 5 together with(9) to obtain the predictions of the activation energy ofa pair of noninteracting kinks at zero applied stress .These are 2 H k = 1 .
80 eV if the dislocation is consid- ered to move along a straight path and 2 H k = 2 .
08 eVfor the curved path of the dislocation [black path inFig. 4(b)]. The latter is in excellent agreement with thevalue 2 H k = 2 .
06 eV obtained in Ref. 13 from the exper-imental data, where the thermal component of the flowstress vanishes. The curvature of the dislocation pathwayand the associated flat maximum of the Peierls barrierwill also give rise to larger curvature of the stress depen-dence of the activation enthalpy and thus to a steepincrease of the flow stress with decreasing temperature,as observed universally in all BCC metals.
V. CONCLUSIONS
We have developed a numerical scheme that providesthe position of the center of the dislocation only fromthe knowledge of the displacements of atoms between therelaxed configuration and the perfect lattice. This repre-sents a generalization of the model devised originally byPeierls and Nabarro . Our model goes beyond thesedevelopments in that it provides two components of thedislocation position, one of which lies in the glide plane,whereas the other is perpendicular to this plane.The position of the dislocation in each of the three pos-sible { } slip planes is identified as the point, wherethe displacement of atoms parallel to the slip directionis half-way between its minima and maxima. This cal-culation is made separately for the chain of atoms aboveand below the lattice site with the dislocation. For eachslip plane, one thus obtains two points that define a linewith possible positions of the dislocation. The pairwiseintersections of the three lines thus obtained define thecorners of a triangle that approximates the position of thedislocation. We have shown that one edge of this trian-gle is ill-defined when the dislocation is very close to theboundary between the adjacent lattice site. However, thedislocation path obtained from the intersections of thetwo remaining edges is smooth and the dislocation po-sitions agree with differential displacement maps as wellas with DFT calculations made by Itakura et al. .We have demonstrated this procedure by calculatingthe path of the dislocation between two minimum en-ergy lattice sites, using the snapshots of atomic positionsobtained from our recent NEB+r calculations made inRef. 17. These calculations have been made for all three { } planes on which the dislocation can move and forzero, positive and negative shear stresses perpendicularto the slip direction. The curvature of the dislocationpath is shown to affect the shape of the Peierls barrier,which becomes steeper close to the beginning and end of the transition path and, at the same time, more flatclose to its maximum. This shape will affect macroscopicpredictions, such as the temperature dependence of theyield stress. In particular, we expect that the activa-tion enthalpy will be a stronger function of the appliedstress, which will give rise to a steeper increase of theflow stresses with decreasing temperature. Acknowledgments
The author thanks Vaclav Vitek for fruitful discussionson the topic. This research was supported by the Marie-Curie International Reintegration Grant No. 247705“MesoPhysDef” and by the Academy of Sciences of theCzech Republic, Project no. RVO:68081723. This workhas been carried out at the Central European Institute ofTechnology (CEITEC) with research infrastructure sup-ported by the project CZ.1.05/1.1.00/02.0068 financedfrom the EU Structural Funds. ∗ Electronic address: [email protected] V. Vitek, R. C. Perrin, and D. K. Bowen, Philos. Mag. ,1049 (1970). W. Xu and J. A. Moriarty, Comp. Mater. Sci. , 348(1998). C. Woodward and S. I. Rao, Philos. Mag. A , 1305(2001). S. L. Frederiksen and K. W. Jacobsen, Philos. Mag. ,365 (2003). J. Li, C.-Z. Wang, J.-P. Chang, W. Cai, V. V. Bulatov,K.-M. Ho, and S. Yip, Phys. Rev. B , 1 (2004). M. Cawkwell, D. Nguyen-Manh, C. Woodward, D. G. Pet-tifor, and V. Vitek, Science , 1059 (2005). L. Ventelon and F. Willaime, J. Comp.-Aid. Mat. Design , 85 (2007). L. Pizzagalli, P. Beauchamp, and H. J´onsson, Philos. Mag. , 91 (2008). R. Gr¨oger, A. G. Bailey, and V. Vitek, Acta Mater. ,5401 (2008). D. Rodney and L. Proville, Phys. Rev. B , 094108(2009). M. Mrovec, D. Nguyen-Manh, C. Els¨asser, and P. Gumb-sch, Phys. Rev. Lett. , 246402 (2011). J. E. Dorn and S. Rajnak, Trans. AIME , 1052 (1964). R. Gr¨oger and V. Vitek, Acta Mater. , 5426 (2008). H. J´onsson, G. Mills, and K. W. Jacobsen,
Classicaland Quantum Dynamics in Condensed Phase Simulations (World Scientific, 1998), chap. 16. Nudged elastic bandmethod for finding minimum energy paths of transitions,pp. 385–404. G. Henkelman and H. J´onsson, J. Chem. Phys. , 9978(2000). R. Gr¨oger and V. Vitek, Model. Simul. Mater. Sci. Eng. , 035019 (2012). R. Gr¨oger and V. Vitek, Acta Mater. , 6362 (2013). R. Peierls, Proc. Phys. Soc. , 34 (1940). F. R. N. Nabarro, Proc. Phys. Soc. , 256 (1947). L. Proville, L. Ventelon, and D. Rodney, Phys. Rev. B ,144106 (2013). L. Ventelon, F. Willaime, E. Clouet, and D. Rodney, ActaMater. , 3973 (2013). J. P. Hirth and J. Lothe,
Theory of dislocations (J.Wiley& Sons, New York, 1982), 2nd ed. M. Itakura, M. Kaburaki, and M. Yamaguchi, Acta Mater. , 3698 (2012). R. Heinrich and W. Schellenberger, Phys. Stat. Sol. B ,81 (1971). L. Dezerald, L. Ventelon, E. Clouet, C. Denoual, D. Rod-ney, and F. Willaime, Phys. Rev. B , 024104 (2014). T. Suzuki, H. Koizumi, and H. O. K. Kirchner, Acta Met-all. Mater. , 2177 (1995). R. Gr¨oger, Ph.D. thesis, University of Pennsylvania (2007),URL http://arxiv.org/abs/0707.3577 . M. Mrovec, R. Gr¨oger, A. G. Bailey, D. Nguyen-Manh,C. Els¨asser, and V. Vitek, Phys. Rev. B , 104119 (2007). G. J. Ackland and R. Thetford, Philos. Mag. A , 15(1987). D. Caillard and J. L. Martin,
Thermally activated mecha-nisms in crystal plasticity (Pergamon Press, 2003). We implicitly assume that all other DOF are relaxed uponspecifying the dislocation position. The minimum of the Peierls barrier V ( ξ ; Σ nonglide ) for agiven applied stress tensor Σ nonglide is shifted to zero. The planes above and below the lattice site with the dislo-cation refer to the two planes adjacent to this site. It is notimportant which of the two planes is considered as being“above” and which “below” this lattice site. Since E ≫ V , the energy in (9) can be closely approxi-mated as 2 H k = √ E R ξ max p V ( ξ ) d ξξ