Free energy generalization of the Peierls potential in iron

Free energy generalization of the Peierls potential in iron

M. R. Gilbert EURATOM/CCFE, Culham Science Centre, Abingdon, Oxfordshire, OX14 3DB, UK    P. Schuck Oak Ridge National Laboratory, Oak Ridge, TN 37831    B. Sadigh Lawrence Livermore National Laboratory, Livermore, CA 94551    J. Marian Lawrence Livermore National Laboratory, Livermore, CA 94551

In body-centered cubic (bcc) crystals, screw dislocations exhibit high intrinsic lattice friction as a consequence of their non-planar core structure, which results in a periodic energy landscape known as the Peierls potential, . The main features determining plastic flow, including its stress and temperature dependences, can be derived directly from this potential, hence its importance. In this Letter, we use thermodynamic integration to provide a full thermodynamic extension of for bcc Fe. We compute the Peierls free energy path as a function of stress and temperature and show that the critical stress vanishes at 700K, supplying the qualitative elements that explain plastic behavior in the athermal limit.

thanks: Corresponding author

Dislocations are ubiquitous line defects that mediate plastic deformation in crystalline materials. In body-centered-cubic (bcc) metals, plasticity is governed by the motion of screw dislocations on close-packed planes. Generally, this motion is understood to occur over a periodic energy landscape known as the Peierls potential . Theoretical descriptions of this potential show that it is very stiff in bcc Fe, leading in some cases to critical stresses (those at which the lattice resistance is suppressed) in excess of one GPa. However, experimentally it is found that the flow stress —the macroscopic equivalent of the critical stress— is roughly one third lower than calculated values. The most convincing explanation for this discrepancy that we possess currently is the contribution of zero-point motion to the Peierls potential at temperatures where quantum effects cannot be neglected proville2012 ().

At low stresses, one can safely assume that the Peierls potential remains unchanged and that slip proceeds via the thermally-activated nucleation of wiggles on the dislocation line, known as kink pairs, and their subsequent sideward relaxation. However, at stresses approaching the critical stress, referred to as Peierls stress at 0K, it is insufficient to consider only the zero stress internal energy to represent the Peierls trajectory. This trajectory is defined as the rectilinear path, denoted by the reaction coordinate , between two equivalent equilibrium states (known as ‘easy core’) on the Peierls potential, which has periodicity , where is the lattice constant. Rodney and Proville proville2009 () showed that, at moderate to high stresses, the core undergoes internal transformations that modify the Peierls energy landscape. This , where is the shear stress resolved on a {110} plane –applied using the Parrinello-Rahman method parinellorahman1981 ()–, is provided in Fig. 1 using the nudged-elastic-band (NEB) method jonssonetal1998 () and a standard semi-empirical interatomic force field developed by Mendelev et al. mendelevhan2003 (). Although this force field yields the non-degenerate screw dislocation core structure, in accordance with density functional theory (DFT) calculations, it also predicts a metastable split core configuration known now to be an artifact ventelon2007 (). The Peierls potential is then defined by the system enthalpy , where is equal to the plastic work associated with the shear strain in the simulation box. For an accurate calculation of , the reaction coordinate must be expressed in terms of the dislocation core position, which is calculated here by matching the atomic displacement field to the Volterra solution for a screw dislocation. For the model used in Fig. 1, is approximately 1250 MPa as first computed by Chaussidon et al. chaussidonetal2006 () (DFT calculations: to MPa mitsu (); ventelon2013 ()).

Figure 1: Peierls potential as a function of stress for the Mendelev interatomic potential for Fe mendelevhan2003 (). represents the (non dimensional) dislocation core position. Note that DFT calculations predict a sinusoidal profile with an amplitude of meV/ ventelon2013 (). The inset to the figure shows the transition path on the (111) plane taken by the dislocation at several values of shown in the main figure.

The inset to Fig. 1 shows the two-dimensional representation of the NEB trajectory at several stresses on the plane using a differential displacement map. The figure shows that, as the shear stress increases, the path approaches one of the atomic rows associated with split core configuration. At zero stress, the transition path is practically rectilinear –consistent with recent DFT calculations mitsu (); ventelon2013 ()– which reflects the structure of the energy landscape for the Mendelev potential mark ().

Tensile tests place the athermal limit of bcc Fe, i.e. the point at which flow occurs without mechanical aid, at various temperatures between 300 and 400K basinski1960 (); mordike1962 (); conrad1962 (); spitzig1970 (); kuramoto1979 (); brunner1991 (). This limit is thought to establish the extent of validity of the classical kink-pair mechanism. Kink-pair energies have been calculated as a function of stress utilizing atomistic and line tension models, all of which make use of a substrate Peierls potential mitsu (); proville2013 (). However, despite its importance, the effect of temperature on the Peierls potential has not yet been addressed. In this Letter we generalize the Peierls enthalpy to finite temperature conditions by calculating the Gibbs free energy of atomistic Fe systems using a combination of periodic (dipole) and cylindrical configurations containing in excess of atoms, as described in Refs. schuck () and rao1998 (). Dislocation segments in length were considered, where is the Burgers vector’s modulus and nm. The reference configurations used here are those calculated at constant stress using the NEB method.

Adding to the natural variables and , results in the isothermal/isobaric ensemble, whose characteristic state function is the Gibbs free energy:

where is the entropy, defined in our pure and periodic systems solely by vibrational contributions. Our objective is to establish the importance of incorporating temperature effects into models based on the standard picture of the Peierls potential and to compute explicity the athermal limit from atomistic calculations. To obtain the Peierls free energy we equate the Peierls transition path to an activated process described by a general configurational nonlinear many-body reaction coordinate carter1989 (). As a first approximation, we first calculate the harmonic free energy of each configuration along the NEB trajectory as:


where are the eigenfrequencies corresponding to eigenvectors pertaining to each NEB replica . However, as Fig. 2 demonstrates, relatively large mean square displacements can be measured already at 100K along the reaction coordinate, particularly within the initial 15% of the Peierls trajectory. Such anharmonic behavior is a manifestation of a pathology of the Mendelev force field, which results in a transition path at 0K that may not reflect the true finite temperature dynamics of the system. has been calculated for comparison nonetheless, and, following Proville et al. proville2012 (), to account for zero-point motion at low temperatures. We shall discuss this correction in later paragraphs.

Figure 2: Time-averaged atomic mean square displacement along the Peierls trajectory (up to 45% of the reaction coordinate) at zero stress. Areas with large indicate anharmonic behavior. The unstressed is shown as a gray dashed line for reference.

Our method to compute full, anharmonic free energies is based on Kirkwood’s approximation to obtain the potential of mean force kirkwood1935 (). Assuming a Hamiltonian of the type , where and are the generalized momenta and coordinates, respectively, one can write the Gibbs free energy as:


where is the canonical partition function, ( is Boltzmann’s constant), and is a constant that represents the integrated contribution of the kinetic energy and the elastic work. If one now extracts the (NEB) trajectory degree of freedom (DOF) from the -dimensional vector , and separate them in eq. 2 we have:


where and are the internal and free energies for the -DOF system defined by generalized coordinates . The force along the trajectory can be evaluated as:


The quantity can be regarded as the configurational partition function of the -DOF system, and, therefore, eq. 4 can be written as:



which, when inserted into eq. 5, results in:


which is a configurational average over all DOF. In other words, the free energy of the constrained system is obtained by integrating the time-averaged total force along the minimum free energy path.

Using eq. 6, we now calculate free energies for the different trajectory points of the NEB calculations in our atomistic systems. Configurational averages are numerically intensive, and require long simulation times to converge (on the order of several ns in our case). The resulting free energies for the unstressed configurations are shown in Fig. 3 at temperatures ranging from 100 to 600K. at 0, 100 and 200K are also included for comparison. Two features are noteworthy at first glance: (i) the characteristic metastability associated with the split core configuration at 0K is lost following constrained equilibration at finite temperatures; (ii) the free energies suffer a marked decrease from 0 to 100K. We find that vanishes completely by 700K at zero stress. Technically, the current force field for Fe is strictly valid above the Debye temperature and below the Curie transition (K). However, Proville et al. have shown that quantum effects in this context are only important below 40K proville2012 (), and so our results over the 100-to-700K temperature range are within the validity margins of the potential.

As mentioned above, the double-‘hump’ shape of the Mendelev interatomic force field as given by the NEB method is a consequence of the metastability of the split core configuration at 0K. By contrast, DFT calculations show that the energy path displays a simpler sinusoidal profile ventelon2007 (); mitsu (). Gilbert et al. have performed a detailed numerical construction of the two-dimensional energy landscape for the Mendelev potential on the {111} plane mark (). Their analysis revealed no justification for the standard NEB trajectory to be favored at 0K. Thus, at finite temperatures the system escapes the NEB path and samples a broader region of phase space until falling into an alternate dynamic path. This path is the one shown in Figure 3 at different temperatures. This behavior explains the long convergence times and the large anharmonicities captured in Fig. 2 as the system samples alternative phase space trajectories. In any case, at 100K both the harmonic and full free energy curves mimic one another for 00.15, which indicates that within that interval the minimum potential and free energy paths are equivalent. At higher temperatures, the range of agreement gradually decreases, as more and more of the barrier is subject to anharmonic effects.

Figure 3: Peierls free energy path at zero stress as a function of temperature. The harmonic free energies at 100 and 200K are also shown as dashed lines for comparison. At 700K and above the barrier is at or below zero.

Next, we study the variation of the free energy barrier with temperature and stress. When displays a weak or no dependence on the stress, the dependence of Peierls enthalpy on is via the temperature-independent plastic work (per unit length). This is shown to be an accurate approximation for MPa (cf. Fig. 1). In such case one need only consider the temperature dependence of the unstressed Peierls trajectory (given in Fig. 3) and subtract to obtain . From this, is measured at each stress and temperature and each value plotted in Fig. 4. We term this approximation (or if referring to free energy barriers). The figure reveals several interesting features. First, the free energy decreases by more than 50 percent from zero to 100K. Subsequently, it decreases gradually until it vanishes. This latter point furnishes the critical stress , i.e. that at which . Second, the curves at different stresses roughly mimic one another within the envelope of the zero stress results, which is an indication that the dynamic path is stable. Much in the manner of Proville et al. proville2012 (), adding zero-point corrections reduces the value of at low temperatures for all stresses. For clarity, this effect is not shown in the figure but will be taken into account when computing the critical stresses.

Figure 4: Variation of the free energy barrier with temperature and stress. Solid lines correspond to , which is valid up to 400 MPa. Scatter points represent i.e. the free energies for the stress dependent potential obtained from NEB. In the latter case, the free energy is strictly zero for MPa.

To verify the assumption of stress-independence for the Peierls potential at 0K below MPa, we have performed free energy calculations on stressed transition paths at stresses 100 to 400 MPa. The results as a function of temperature are shown as scatter points in Fig. 4. As shown, the agreement between and the stressed configurations () is reasonable up to 300 MPa. From the data points presented in Fig. 4, we can extract the values of at each temperature. These are shown in Fig. 5, where the corresponding harmonic values are displayed for comparison. At 0K, including zero-point motion results in MPa, which we use as a common point (joined by dashed lines) for all the curves.

Figure 5: Variation of the critical stress with temperature from full and harmonic free energy calculations. Solid lines represent classical calculations, whereas the dashed ones symbolize the quantum corrections at low temperature. In the classical limit, at K.

The critical stress vanishes at K, which represents the athermal limit within our model. This is higher than experimentally observed but our results show that the mere consideration of temperature effects leads to important changes in the dynamic picture of screw dislocations. As an example, at 100K, MPa, which represents a 60% reduction with respect to the value of at 0K.

The data provided in Fig. 5 are a central result of this paper and reveal two main behaviors related to the Peierls potential. As referred to earlier, tensile tests in pure Fe show that the temperature dependence of the flow stress is characterized by the critical stress at (very) low temperatures ( K) and the athermal limit at high temperatures. In the absence of quantum corrections, atomistic calculations, even of the most accurate kind, fail to predict the lower temperature limit, while there is no numerical work available providing information for the higher one. The present calculations show that a formal treatment of the free energy may account for dramatic reductions in both limits using conventional interatomic force fields. Indeed, these calculations provide a closed set of data for defining a temperature-dependent substrate potential to be used in higher level models such as line tension, line-on-substrate, kinetic Monte Carlo, etc. We believe the implications of this work to be of importance to all bcc metals. It is worth emphasizing that the non-sinusoidal nature of from the force field employed here is not a weakening aspect of this work because finite-temperature trajectories sample paths in phase space that are not affected by the existence of the split core configuration at 0K. This may explain why most molecular dynamics simulations of screw dislocation motion using the present interatomic potential show only correlated (in the sense of Gordon et al. gordon2010 ()) formation of kink pairs domain2005 (); gilbert2011 () at finite temperatures.

To conclude, we have presented a free energy map of the screw dislocation core transition on planes. The calculations have been done using constrained reaction coordinate dynamics and reveal a drastic reduction in free energy barrier and critical stress with increasing temperature. Our results can serve as yet another platform from which to interpret the discrepancies observed between atomistic simulations and macroscopic flow stress measurements.

The authors gratefully acknowledge helpful discussions with D. Rodney, S. Dudarev, and P. Derlet. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. This work was part-funded by the RCUK Energy Programme under grant EP/I501045 and the European Communities under the contract of Association between EURATOM and CCFE. The views and opinions expressed herein do not necessarily reflect those of the European Commission.


  • (1)
  • (2) L. Proville, D. Rodney, M.-C. Marinica, Nat. Mater. 11, 845 (2012).
  • (3) D. Rodney and L. Proville, Phys. Rev. B 79, 094108 (2009).
  • (4) M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
  • (5) H. Jonsson, 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.
  • (6) M. I. Mendelev, S. Han, D. J. Srolovitz, G. J. Ackland, D. Y. Sun, and M. Asta, Philos. Mag. 83, 3977 (2003).
  • (7) L. Ventelon and F. Willaime, J. Comput.-Aided Mater. Des. 14 85 (2007).
  • (8) J. Chaussidon, M. Fivel, D. Rodney, Acta Mater. 54, 3407 (2006).
  • (9) M. Itakura, H. Kaburaki, and M. Yamaguchi, Acta Mater. 60, 3698 (2012).
  • (10) L. Ventelon, F. Willaime, E. Clouet, D. Rodney, Acta Mater. 61, 3973 (2013).
  • (11) M. R. Gilbert, P. M. Derlet and S. L. Dudarev,
  • (12) Z. S. Basinski, and J. W. Christian, Australian Journal of Physics 13, 299 (1960).
  • (13) B. L. Mordike and P. Haasen, Philos. Mag. 7, 459 (1962).
  • (14) H. Conrad, S. Frederick, Acta Metall. 10, 1013 (1962).
  • (15) W. A. Spitzig and A. S. Keh, Acta Metall. 18, 611 (1970).
  • (16) E. Kuramoto, Y. Aono, and K. Kitajima, Philos. Mag. A 39, 717 (1979).
  • (17) D. Brunner and J. Diehl, Phys. Status Solidi A 124, 155 (1991).
  • (18) L. Proville, L. Ventelon, and D. Rodney, Phys. Rev B 87, 144106 (2013).
  • (19) P. C. Schuck, J. Marian, J. B. Adams, and B. Sadigh, Philos. Mag. 89, 2861 (2009).
  • (20) S. Rao, C. Hernandez, J. P. Simmons, T. A. Parthasarathy, and C. Woodward, Philos. Mag. A 77, 231 (1998).
  • (21) E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. Lett. 156, 472 (1989).
  • (22) J. G. Kirkwood, J. Chem. Phys. 3 300 (1935).
  • (23) A. Gordon, T. Neeraj, Y. Li, and J. Li, Modell. Simul. Mater. Sci. Eng. 18, 085008 (2010).
  • (24) C. Domain and G. Monnet, Phys. Rev. Lett. 95, 215506 (2005).
  • (25) M. R. Gilbert, S. Queyreau, and J. Marian, Phys. Rev. B 84 174103 (2011).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description