First Principles Based Multiscale Modeli

First Principles Based Multiscale Modeling of Single Crystal Plasticity: Application to BCC Tantalum - CaltechTHESIS
CaltechTHESIS
A Caltech Library Service
About
Browse
Deposit an Item
Instructions for Students
First Principles Based Multiscale Modeling of Single Crystal Plasticity: Application to BCC Tantalum
Citation
Wang, Guofeng
(2002)
First Principles Based Multiscale Modeling of Single Crystal Plasticity: Application to BCC Tantalum.
Dissertation (Ph.D.), California Institute of Technology.
doi:10.7907/5nyn-ct36.
Abstract
In principle, the macroscopic plasticity properties of crystalline materials are derivable from the physical processes involving dislocations and interactions between dislocations with other defects. However, a quantitative theory of plasticity based on the dislocation mechanism requires crossing multiple length and time scales. To accommodate these requirements, we developed a multiscale approach for modeling crystalline solids. In this thesis, to establish the connections between simulations in different length and time scales, I mainly focus on identifying and determining the importance and influence of various unit processes involving the dislocations through atomic level simulations. These unit processes in turn play a major role in modeling the single crystal plasticity.
Key Results from Atomistic Simulations
Dislocation core structure and core energy: Using the first-principles qEAM force field (FF), we determine the core energy for 1/2a<111> screw dislocation and 1/2a<111> edge dislocation in bcc Ta. We find that the core energy of edge dislocation is 1.77 times higher than that of screw dislocation. This ratio (1.77) is a fundamental material property used as input to the macroscopic model. Furthermore, we find that the central 12 atoms closest to the 1/2a<111> screw dislocation line have distinguishably higher atomistic strain energy than the other atoms. Thus, we arrive at a physical definition of dislocation core.
Screw dislocation mobility: In this thesis, we proposed a new method to investigate dislocation mobility by analyzing the process of migration of a screw dislocation dipole. The new method is based on the energy distribution at the atomistic scale and is used to calculate the Peierls potential barrier and Peierls stress for dislocation continuous motion. The calculated Peierls stress is in good agreement with results obtained using other method. Simulating dislocation motion at finite temperatures (from 20 K to 300 K), we find that the activation energy for dislocation motion is about 6 times lower than computed at 0.001 K. Our results suggest that the decrease in the correlation between neighboring segments in the dislocation line accounts for the decrease of activation energy. We observe that the formation of kink pair along the dislocation line enhances the dislocation mobility. This verifies the traditional belief that the screw dislocation in bcc metals moves by first kink pair nucleation and subsequently lateral movements of kinks along the dislocation.
Kinks in screw dislocations: To bridge the atomistic process of dislocation motion with continuum model, we accurately calculate the material properties, such as kink pair formation energy and effective kink pair length, using atomic level simulations. In detailed structural analysis, we discover the substructures of different kinks when the screw dislocation core is asymmetric. There are only two kinds of elementary kinks in the dislocation and the others are the composite kinks consisting of an elementary kink and one or two flips. Based on these findings, we further explain the observed trend of the formation energy and mobility of different classes of kinks. (Note: Similar trend and conclusion could have been found in earlier studies but not mentioned by the authors of those papers.)
In summary, we have used quantum mechanics based interaction potentials to investigate the unit processes that play important role in single crystal plasticity and verified the findings using the quantitative results obtained from the atomic level simulation in a macroscopic model for single crystal plasticity.
Item Type:
Thesis (Dissertation (Ph.D.))
Subject Keywords:
Materials Science; Computer Science
Degree Grantor:
California Institute of Technology
Division:
Engineering and Applied Science
Major Option:
Materials Science
Minor Option:
Computer Science
Thesis Availability:
Public (worldwide access)
Research Advisor(s):
Goddard, William A., III
Thesis Committee:
Goddard, William A., III (chair)
Cagin, Tahir
Fultz, Brent T.
Haile, Sossina M.
Johnson, William Lewis
Wang, Zhen-Gang
Defense Date:
30 April 2002
Record Number:
CaltechTHESIS:11132009-112545862
Persistent URL:
DOI:
10.7907/5nyn-ct36
ORCID:
Author
ORCID
Wang, Guofeng
0000-0001-8249-4101
Default Usage Policy:
No commercial reproduction, distribution, display or performance rights in this work are provided.
ID Code:
5374
Collection:
CaltechTHESIS
Deposited By:
Tony Diaz
Deposited On:
17 Nov 2009 21:56
Last Modified:
23 Aug 2022 22:49
Thesis Files
PDF
- Final Version
See Usage Policy.
5MB
Repository Staff Only:
item control page
CaltechTHESIS is powered by
EPrints 3.3
which is developed by the
School of Electronics and Computer Science
at the University of Southampton.
More information and software credits
First Principles Based Multiscale Modeling of
Single Crystal Plasticity:
Application to BCC Tantalum

Thesis by

Guofeng Wang

In Partial Fulfillment of the Requirements
for the Degree of
Doctor of Philosophy

California Institute of Technology
Pasadena, California

2002
(Defended April 30, 2002)

Guofeng Wang

iii

Acknowledgement
I would like to thank all the people who helped me in my study, research, and life
at Caltech. I appreciate their kindness and assistance. Especially, I would like to thank the
following people for their contributions to this thesis.

First, I would like to thank my advisor, Prof. William A. Goddard III, for his
support, encouragement and guidance in the last five years. Bill is a great scientist, who is
knowledgeable, creative, and enthusiastic. He patiently educated and tutored me in
scientific thinking and English writing and speaking. He sets a fine example for me to
follow in my future career.

Second, I would like to thank Dr. Tahir <;agm and Dr. Alejandro Strachan. They
are my project directors and good friends. We discussed research works deeply and
thoroughly and I learned many techniques from them. Furthermore, they also gave me
many suggestions on my life and career. It is an honor for me to acknowledge Prof.
William L. Johnson, Prof. Brent Fultz, Prof. Zhen-Gang Wang and Prof. Sossina Haile
for their attendance at both my candidacy exam and thesis defense. They are the best
committee members. They gave me much advice on my work and personal development.
During my stay at Caltech, I benefited a lot from collaboration with other current and
former group members. I would like to thank everyone at MSC. In particularly, Ryan
Martin, Yue Qi, Darryl Willick, Rick Muller, Peter Schultz (Sandia National Laboratory),
Hyon-Jee Lee and Darren Segall (MIT).

Last but not least, I dedicate this thesis to my parents, my wife and my daughter. I
would like to express my gratitude to my parents, who have been supporting my
education both financially and emotionally. My wife, Fang Li, gives me love, support and
encouragement. My daughter Eilleen Wang, who is definitely an angel, brings us
happiness and makes us proud.

Abstract
In principle, the macroscopic plasticity properties of crystalline materials are
derivable from the physical processes involving dislocations and interactions between
dislocations with other defects. However, a quantitative theory of plasticity based on the
dislocation mechanism requires crossing multiple length and time scales. To
accommodate these requirements, we developed a multiscale approach for modeling
crystalline solids. In this thesis, to establish the connections between simulations in
different length and time scales, I mainly focus on identifying and determining the
importance and influence of various unit processes involving the dislocations through
atomic level simulations. These unit processes in tum play a major role in modeling the
single crystal plasticity.

Key Results from Atomistic Simulations
Dislocation core structure and core energy: Using the first-principles qEAM force field (FF),
we determine the core energy for 1/2a screw dislocation and 1/2a edge
dislocation in bee Ta. We find that the core energy of edge dislocation is 1.77 times higher
than that of screw dislocation. This ratio (1.77) is a fundamental material property used as
input to the macroscopic model. Furthermore, we find that the central 12 atoms closest to the
l/2a screw dislocation line have distinguishably higher atomistic strain energy than the
other atoms. Thus, we arrive at a physical definition of dislocation core.
Screw dislocation mobility: In this thesis, we proposed a new method to investigate
dislocation mobility by analyzing the process of migration of a screw dislocation dipole. The
new method is based on the energy distribution at the atomistic scale and is used to calculate

the Peierls potential barrier and Peierls stress for dislocation continuous motion. The
calculated Peierls stress is in good agreement with results obtained using other method.
Simulating dislocation motion at finite temperatures (from 20 K to 300 K), we find that the
activation energy for dislocation motion is about 6 times lower than computed at 0.001 K.
Our results suggest that the decrease in the correlation between neighboring segments in the
dislocation line accounts for the decrease of activation energy. We observe that the formation
of kink pair along the dislocation line enhances the dislocation mobility. This verifies the
traditional belief that the screw dislocation in bee metals moves by first kink pair nucleation
and subsequently lateral movements of kinks along the dislocation.
Kinks in screw dislocations: To bridge the atomistic process of dislocation motion with
continuum model, we accurately calculate the material properties, such as kink pair
formation energy and effective kink pair length, using atomic level simulations. In detailed
structural analysis, we discover the substructures of different kinks when the screw
dislocation core is asymmetric. There are only two kinds of elementary kinks in the
dislocation and the others are the composite kinks consisting of an elementary kink and_one
or two flips. Based on these findings, we further explain the observed trend of the formation
energy and mobility of different classes of kinks. (Note: Similar trend and conclusion could
have been found in earlier studies but not mentioned by the authors of those papers.)

In summary, we have used quantum mechanics based interaction potentials to
investigate the unit processes that play important role in single crystal plasticity and
verified the findings using the quantitative results obtained from the atomic level
simulation in a macroscopic model for single crystal plasticity.

vii

Table of Contents

Acknowledgement ..................................................................................iii
Abstract ................................................................................................ v
Table of Contents ................................................................................. vii

Chapter 1 Introduction

1.1 Dislocations and plasticity.......................................................... 1
1.2 Features in plasticity of bee metals ................................................ 4
1.3 Computer modeling of material plasticity ........................................ 6
1.4 References .............................................................................. 11

Chapter 2 Atomistic simulation methods
2.1 Overview ............................................................................... 14

2.2 Molecular dynamics (MD) simulations .......................................... 15
2.3 Embedded atoms model (EAM) force field ..................................... 19
2.3.1 Physical foundation ....................................................... 19
2.3.2 The qEAM FF ............................................................. 22

2.3.3 Validation of the qEAM FF .•........................................... 24
2.4 References............................................................................ 32

Chapter 3 Core structure and core energy of 1/2a screw dislocation in Ta
3.1 Overview ............................................................................... 37

3.2 Dislocation core structure ..........................................................38

Vlll

3.2.1 Construction of the dislocation quadrupole .........................38
3.2.2 The differential displacement (DD) map ............................. 40
3.2.3 Polarization of the dislocation ......................................... .41
3.2.4 Comparison to other calculations ..................................... .45
3.3 Dislocation core energy .............................................................. 47
3.3.1 Atomistic approach ....................................................... 47
3.3.2 Continuum approach ..................................................... 51

3.4 Dislocation core structure revisited ...............................................57
3.4.1 Dislocation core energy variations with its polarization ..........57
3.4.2 Force field re-parameterization ........................................61
3.4.3 Generalized stacking-fault energy (y) surface .......................64
3.5 Conclusion ............................................................................. 68
3.6 References .............................................................................. 70

Chapter 4 Peierls energy barrier and Peierls stress of 1/2a screw dislocation in
Ta
4.1 Overview ............................................................................... 72

4.2 Dislocation dipole migration and annihilation process .......................73
4.3 Peierls energy barrier and Peierls stress from dislocation dipole
migration ............................................................................... 80

4.4 Twinning/Anti-twinning asymmetry ............................................. 84
4.5 Dislocation motion at finite temperatures .......................................87
4.6 Dislocation motion by nucleating kinks .......................................... 96
4. 7 Conclusion ............................................................................. 99

4.8 Appendix: Periodic boundary for the simulation cell containing a screw
dislocation dipole .................................................................... 101
4.9 References ............................................................................ 106

Chapter 5 Flips and kinks on 1/2a screw dislocation in Ta
5.1 Overview .............................................................................. 107
5.2 Simulation model ................................................................... 109

5.2.1 Construction of simulation model .................................... 109
5.2.2 Boundary conditions of simulation models ........................ 111
5.3 Multiplicity of flips and kinks .................................................... 114
5.3.1 Equilibrium dislocation core structures ............................ 114
5.3.2 Flips ........................................................................ 118
5.3.3 Isolated kinks ............................................................. 118

5.4 Kink formation energy calculations ............................................. 120
5.4.1 Formation energy of the isolated kink .............................. 120
5.4.2 Formation energy of kink pairs ....................................... 124
5.5 Kink migration energy calculations ............................................. 127
5.5.1 Kink migration energy ................................................. 127

5.5.2 Relative mobility of kinks .............................................. 129
5.6 Structural analyses .................................................................. 134
5.6.1 Overview .................................................................. 134

5.6.2 Structural analysis of flips ............................................. 135
5.6.3 Structural analysis of kinks ........................................... 146
5.6.4 Determination of geometrical parameters .......................... 160

5. 7 Conclusion ........................................................................... 164
5.8 References ............................................................................ 166

Chapter 6 A multiscale approach for modeling crystalline solids
6.1 Overview .............................................................................. 169
6.2 Unit processes ........................................................................ 171

6.2.1 Dislocation mobility: double kink formation and thermally
activated motion of kinks .............................................. 172
6.2.2 Dislocation interactions: obstacle-pair strength and obstacle
strength .................................................................... 176

6.2.3 Dislocation evolution: multiplication and attrition ............... 182
6.3 Atomistic modeling of dislocation properties ................................. 184
6.3.1 Core energy of 1/2a screw and edge dislocations .......... 186
6.3.2 Kink pair energy and nucleation length ............................ 193
6.4 Experiment, validation and predication ....................................... 198
6.5 Conclusion ........................................................................... 205
6.6 Comment .............................................................................. 206
6. 7 References ............................................................................ 207

Chapter 1

Chapter 1 Introduction
1.1 Dislocations and plasticity
Dislocations 1•2 •3 are line defects in the atomic arrangement of a crystalline
material. Since dislocation is a disorder m the crystalline system, its presence m a
material increases the internal energy, electrical conductivity, and hardness, and
influences many other physical properties. Among all influences caused by dislocations,
people are most interested in the role that dislocations play in the plastic deformation of
materials. Many experimental and theoretical studies have established the belief that
dislocations are the primary agents of plasticity, i.e., plastic deformation proceeds by the
generation and movement of dislocations. It is also firmly established that the
macroscopic plasticity properties of crystalline materials are derivable, at least m
principle, from the behavior of dislocations and their interactions with other defects 4 .

The following provides a brief introduction to dislocations. For a further reading,
please consult the references (Refs. 5, 6 and 7).

(1) Dislocation type. There are two basic types of dislocation in the crystalline
materials: edge dislocation and screw dislocation. Figure 1-1 shows the descriptive model
of the basic geometry of an edge (Figure 1-1 (a)) and a screw (Figure 1-1 (b)) dislocation.
In an edge dislocation, an extra plane of atoms is inserted in the crystal but not extending
through all of the crystal by ending in the dislocation line as illustrated in Figure 1-1 (a).
Figure 1-1 (b) shows a screw dislocation originates from a shift of one atom in the lattice
with respect to a perfect arrangement and can be described as a single surface helicoid,

Chapter 1

rather like a spiral staircase.

In the most general case, the dislocation (called mixed

dislocation) has a mixed edge and screw character.

FORCE

I-

1-1-

-A I

I-

~~ ,1~

(a)

r,. /

- , , ~ ., ., - , , -- ·vV ./ ., .,, ,, ., ✓ .,,

I-

.,,.6'

.FI

-•

,I'

.,~

1)/f./_/

,...

., .,

",,,,l'

,, , , , •

~,,,,.

.., ,,

(;
L,,

I-'

.....
...

... ... ~

....~~ ~

. ;11

~,,J>J>.#1-i.,,-J'

I-

,., I-

,.,

~I-

- - - - , ., , ,, v",,
--- I~.,

.....

FORCE

(b)

Figure 1-1. Descriptive models of the edge dislocation and screw dislocation. (a) In an
edge dislocation, an extra plane of atoms is present in the lattice structure of the crystal.
In this case the extra plane is found adjacent to area A. (b) In a screw dislocation, a plane
of atoms forms a step in the crystal surface. Other atoms can then line up against this step
as the crystal grows.

(2) Burgers vector. Burgers vector b is the fundamental quantity defining an
arbitrary dislocation. Its atomistic definition follows from a Burgers circuit around the
dislocations in the real crystal, which is illustrated in Figure 1-2 for an edge dislocation.
In Figure 1-2 (a), if making a closed circuit from lattice point to lattice point (or atom to
atom) that encloses the dislocation, we obtain a closed chain of the base vectors defining
the lattice. However, making exactly the same chain of base vectors in a perfect reference

Chapter 1

lattice (Figure 1-2 (b)), we would obtain a chain not closed. The vector needed for
closing the circuit in the reference crystal is the Burgers vector b.

i.
.i

(a)

(b)

Figure 1-2. (a) Burgers circuit around an edge dislocation, (b) the same circuit in a
perfect reference crystal; the closure failure is the Burgers vector.

The following are two important rules for dislocation Burgers vector.
(a) The Burgers vector of an edge dislocation is normal to the line of the dislocation.

(b) The Burgers vector of a screw dislocation is parallel to the line of the dislocation.

For a mixed dislocation, the dislocation line may lie at an arbitrary angle to its
Burgers vector. However, the Burgers vector of the dislocation is always the same and
independent of the position of the dislocation.

(3) Slip and the Schmid law.

There are two basic types of dislocation movement: glide (or conservative)
motion, in which the dislocation moves in the surface which contains both its line and

Chapter I

Burgers vector; and climb or (non-conservative) motion, in which the dislocation moves
out of the glide surface normal to the Burgers vector. The concept of slip, providing a
valuable understanding of the structure of the dislocation, is the most important
manifestation of glide.

Plastic deformation in a crystal occurs by the sliding or successive displacement
of one plane of atoms over another on the slip planes. Discrete blocks of crystal between
two slip planes remain undistorted during the slip. Further deformation occurs either by
more movement on existing slip planes or by the formation of new slip planes. The slip
plane for a dislocation is normally the plane with the highest density of atoms and the
direction in the slip is the direction of the slip plane in which atoms are most closely
spaced. A slip plane and a slip direction in the plane constitute a slip system.

A characteristic shear stress is required for dislocation to slip. Arduous
experiments on the relative orientation between the required shear stress and slip system
for a dislocation lead to the famous Schmid law of the critical resolved shear stress
(CRSS)8. The Schmid law states that the dislocation can slip in a slip system when the
shear stress, resolved on the slip plane and in the slip direction, reaches the critical
resolved shear stress (CRSS).

1.2 Features in plasticity of bee metals

For body-centered cubic metals (e.g., iron, molybdenum, tantalum, vanadium,
chromium, tungsten, niobium, sodium, and potassium), there are certain macroscopic

Chapter 1

features common to the low temperature deformation behavior that distinguish the whole

group from fee and hep metals and alloys .

These features include

a rapid increase of the yield and flow stresses with decreasing temperature,

a marked sensitivity of the stress to the imposed strain rate,

a rather small and not very temperature sensitive work-hardening rate,

a sensitivity to small amounts of impurity or solute, particularly, interstitial
solutes,

a tendency in many cases to brittle cleavage fracture at low temperatures,

a complete breakdown of the Schmid law of critical resolved shear stress.

No doubt, associated to these features are the distinguished microscopic physical
processes of dislocations in the bee metals. The in situ high-voltage electron microscope
study at low temperature finds for bee metals the plastic deformation is characterized by
the slow movement of long screw dislocations and fast movement of mixed
dislocations 10 . This leads to the assertion that the mobility of screw dislocations governs
the low temperature deformation behavior for bee metals. Furthermore, the kink pair
mechanism 11 , assuming the screw dislocation in bee metals moves by the kink pair
nucleation and subsequently lateral motion of the component kinks, can successfully
account for the rapid increase of flow stress with decreasing temperature 12- 14 . The most
interesting feature in bee metal is the asymmetry of the slip, which contradicts the
Schmid law. In bee metals (for instances, iron and silicon-iron alloys 15 ,1 6 , tungsten 17 •18
niobium 19 ·20 , tantalum 18 •21 , and molybdenum 21 ' 22 ), the shear stress to move a dislocation

Chapter 1

lying in a slip plane in one direction is not the same as the shear stress required to move it
in the opposite direction in the same plane.

To understand these nontrivial features of plastic deformation in bee metals at low
temperatures, we performed the accurate and systematic simulations for the l/2a
screw dislocation in bee Ta single crystal and summarized the results in this thesis.

1.3 Computer modeling of material plasticity
Following the postulations of dislocations in 1934 1-3 , there have been several
waves of activity in dislocation studies (see reviews 23, 24). The isotropic elastic field
theories of dislocations and interactions among them were developed in the 1940s.
Anisotropic elastic theory, pileup theory, direct observations of dislocations in
transmission electron microscopy and work hardening theory were developed in the
1950s and 1960s. Extended dislocation arrays and the advent of atomistic computer
simulations appeared in the 1960s and 1970s. And, since then, there has been the
refinement in the details of dislocation interactions and core structures, extensive work on
thin films, and computer simulations at several size scales.

With the fast development of the computer power and algorithms, the structures
and behaviors of dislocations could be simulated using the physics and chemistry realistic
models to attain a fundamental understanding of the elementary process of dislocation
slip 25 -30 . However, it is still prohibitive for us to quantitatively derive the macroscopic
material plasticity based on the microscopic dislocation mechanism, because of the need
to trace the evolution of a large number of interacting dislocations over long periods of

Chapter 1

time. On the other hand, many equations of crystal plasticity used for continuum
modeling have been developed to handle the multiplicity and complexity of describing
the mechanisms of dislocation motion and interactions. However, most of these current
continuum equations are phenomenological and largely disconnected from the physics of
the underlying dislocation behavior.

To bridge the existing gap between dislocation physics and crystal plasticity, the
strongly connected multiscale simulations, which are over multiple size scales and time
scales, are necessary31 -39 . In an embedded (or hybrid) multiscale simulation, different
regions are treated in different ways. The region of the greatest interest is simulated using
quantum mechanics (QM) or molecular dynamics (MD), while the region, in which the
atoms move collectively, can be simulated using finite element method (FEM) or field
theory. Another model for the multiscale simulation is the hierarchical informed model.
In this model, separate simulations are carried out at different size scales ranging from

Angstrom in QM regime to meters in continuum materials. At each size scale, some
important physical parameters are extracted from detailed simulations and are input into
the next level simulation with a larger size and time scale.

To better understand the concept of multiscale simulation, Figure 1-3 shows that
the size and time scale of the physical processes could be studied by different simulation
methods. Generally and roughly, the simulation methods are classified into four regimes.
(1) Quantum Mechanics (QM) simulation.

In quantum mechanics calculations, we regularly need to solve the Schrodinger
equation for one or more particles (in most cases, electrons) to obtain the energy and

Chapter 1

force of the system. The popular methods in this field are the Hartree-Fock (HF)
method40 ,4 1 and density functional theory (DFT) method42 ' 43 . The physical system can be
efficiently treated with QM calculations are within the size no larger than tens of
nanometers and over the period no longer than nanoseconds.

Time

I Electrons

Atoms

Grains

Minutes
Seconds

Microsec
1 IJ'

QM '

Picosec

..,
,J

MESO'

11'

Continuum

femtosec

+ ..,

Hours

Nanosec

Material

.4~

..,

,u0
....... Distance

micro

meters

Figure 1-3. Model for the hierarchical informed multiscale simulation.

(2) Molecular dynamics (MD) simulation.

Systems of many interacting atoms or molecules can be studied classically by
solving Newton's equations of motion in an MD simulation. The MD simulation consists
essentially of integrating the equations of motion of the system numerically. Therefore it
simulates the system as it develops over a period of time. In the simulation, the system
moves in the phase space along its physical trajectory as determined by the equations of
motion. Currently, a highly paralleled MD simulation can handle a system over millions

Chapter I

of atoms (or molecules) and over a period of time of microseconds. In Chapter 2, we will
discuss MD simulation methods further.

(3) MESO scale simulation.

To study the subjects such as grain growth or dislocation pattern, the mesoscale
simulation techniques ignoring the atomistic details of the system are desired. There are
many ways to reduce the complexity of the system of interest. The Kinetic Monte Carlo
(KMC) method44 is widely used to simulate a large system over a rather long time based
on the known mechanisms.
(4) Continuum simulation.

In this regime, the material is considered as a continuum media. The fully
developed theories, such as statistical mechanics, kinetic mechanics, and continuum
mechanics, are employed to investigate the material properties in much larger spatial and
temporal scales (for example, Ref. 45).

In this study, we adopt the hierarchical informed model to simulate the single
crystal plasticity for Ta. The red lines in Figure 1-3 show our approach to cross over the
simulations of electrons (QM), atoms (MD), grains (MESO), and material (Continuum).
Our multiscale approach for modeling Ta crystalline solids consists of three hierarchical
parts.

(1) Derive the atomistic interaction potential for Ta based on the data obtained

from the accurate quantum mechanics calculation,

Chapter 1

(2) Predict the properties and behaviors of dislocations m the atomistic
simulations using the derived first-principles potential,

(3) Describe the material plasticity in the kink pair mechanism based mesoscopic
model with the input of the predicted atomistic-level dislocation properties.

This thesis (from Chapter 3 to Chapter 5) will focus on the work of simulating,
identifying and predicting the physical processes of dislocations in atomistic level
simulations, i.e., the part (2) of the whole multiscale simulation approach. The part (1) of
the approach will be briefly described in Chapter 2. Chapter 6 reports the predicted
results of the developed approach exercised to describe the mechanical response of highpurity Tantalum single crystals.

Chapter 1

1.4 References

1. G. I. Taylor, Proc. R. Soc. Lond. A, 145, 362 (1934).
2. E. Orowan, Z. Phys., 89, 605 (1934).
3. M. Polanyi, Z. Phys., 89, 660 (1934).
4. V. V. Bulatov, N. Tang and H. M. Zbib, MRS Bull., 26(3), 191 (2001).
5. D. Hull, Introduction to Dislocations, 2nd Ed., (Pergamon, Oxford, 1975).
6. J. P. Hirth and J. Lothe, Theory of dislocations (Wiley, New York, 1982).
7. Series of Dislocation in Solids, Edited by F. R. N. Nabarro, (North Holland
Publishing Company, Armsterdam).
8. E. Schmid, Proc. Intl. Congr. Appl. Mech., Delft, 342 (1924).
9. J. W. Christian, Metall. Trans. A, 14A, 1237 (1983).
10. F. Louchet and L. P. Kubin, Philos. Mag. A, 39, 433 (1979).
11. A. Seeger and P. Schiller, Physical Acoustics, edited by W. P. Mason (Academics,
New York 1966), Vol. 3A, p. 361.
12. M. S. Duesbery, Acta Metall., 31, 1747 (1983).
13. M. S. Duesbery, Acta Metall., 31, 1759 (1983).
14. M. S. Duesbery and Z. S. Basinski, Acta Metall. Mater., 41, 643 (1993).
15. T. Taoka, S. Takeuchi, and E. Furubayashi, J. Phys. Soc. Japan, 19, 701 (1964).
16. S. Sestak and N. Zarubova, Phys. Stat. Sol., 10, 239 (1965).
17. A. S. Argon and S. R. Maloof, Acta Met., 14, 1449 (1966).
18. D. Hull, J. F. Byron, and F. W. Noble, Can. J. Phys., 45, 1091 (1967).
19. D. K. Bowen, J. W. Christian, and G. Taylor, Can. J. Phys., 45, 903 (1967).

Chapter 1

20. R. A. Foxall, M. S. Duesbery, and P. B. Hirsch, Can. J. Phys., 45,607 (1967).
21. P. J. Sherwood, F. Guiu, H. C. Kim, and P. L. Pratt, Can. J. Phys., 45, 1075 (1967).
22. D. F. Stein, Can. J. Phys., 45, 1063 (1967).
23. F. R. N. Nabarro, Theory of Crystal Dislocations, (Oxford University Press, Oxford,
1967).
24. J.P. Hirth, Metall. Trans., 16A, 2085 (1985).
25. V. V. Bulatov, S. Yip and A. S. Argon, Philos. Mag. A, 72,453 (1995).
26. V. Vitek, Cryst. Lattice Defects, 5, 1 (1974).
27. W. Xu and J. A. Moriarty, Phys. Rev. B, 54, 6941 (1996).
28. L. H. Yang, P. Soderlind, and J. A. Moriarty, Philos. Mag. A, 81, 1355 (2001).
29. S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett., 84, 1499 (2000).
30. C. Woodward and S. I. Rao, Philos. Mag. A, 81, 1305 (2001).
31. V. V. Bulatov et al., Nature, 391, 669 (1998).
32. J. A. Moriarty, W. Xu, P. Soderlind, J. Belak, L. H. Yang and J. Zhu, J. Eng. Mater.
Tech., 121, 120 (1999).
33. A. M. Cuitino, L. Stainier, G. Wang, A. Strachan, T. c;agm, W. A. Goddard, and M.
Ortiz, A Multiscale Approach for Modeling Crystalline Solids, J. Comput.-Aided
Mater. Design, to be published.
34. V. V. Bulatov and L. P. Kubin, Curr. Opin. Solid State Mater. Sci., 3 (6), 558 (1998).
35. R. Phillips, Curr. Opin. Solid State Mater. Sci., 3 (6), 526 (1998).
36. G. H. Campbell, S. M. Foiles, H. C. Huang, D. A. Hughes, W. E. King, D. H. Lassila,
D. J. Nikkel, T. D. de la Rubia, J. Y. Shu, and V. P. Smyshlyaev, Mater. Sci. Eng. A,
251 (1-2), 1 (1998).

Chapter 1

37. R. Phillips, D. Rodney, V. Shenoy, E. Tadmor, and M. Ortiz, Model. Simul. in Mater.
Sci. Eng., 7 (5), 769 (1999).
38. M. I. Baskes, Curr. Opin. Solid State Mater. Sci., 4 (3), 273 (1999).
39. R. Phillips, Crystals, Defects and Microstructures -

Modeling Across Scales,

(Cambridge University Press, 2001).
40. D.R. Hartree, Proc. Camb. Phil. Soc., 24, 89 (1928).
41. F. Fock, Z. Phys., 61, 126 (1930).
42. P. Hohenberg and W. Kohn, Phys. Rev., 136, B864 (1964).
43. W. Kohn and L. J. Sham, Phys. Rev., 140, Al 133 (1965).
44. V. V. Bulatov, J. F. Justo, W. Cai, S. Yip, A. S. Argon, T. Lenosky, M. de Koning,
and T. D. de la Rubia, Philos. Mag. A, 81 (5), 1257 (2001).
45. A. M. Cuitino, L. Stainer and M. Ortiz, Micromechanical modeling of hardening, rate
sensitivity and thermal softening in bee single crystals, Journal of the Mechanics and
Physics of Solids, 2001.

Chapter 2

Chapter 2 Atomistic Simulation Methods
2.1 Overview
Computer simulations have been extensively used in the last several decades and
become an indispensable part of scientific research. Providing detailed linkages between
microscopic and macroscopic properties for the interested system, computer simulations
help us to interpret and design experiments. We could simulate the response of a system
and compare to the experimental observed values to understand the underlying physical
processes. We also could simulate the system under conditions where experiments have
not been performed or cannot be performed easily. Computer simulations in these cases
are able to give the detailed microscopic information that is useful for designing better
experiments. Among various simulation methods, molecular dynamics (MD) is one of the
most widely used simulation methods for studying the properties of liquid, solids and
molecules 1. In MD simulation, the motion of individual particles (atoms or molecules) is
modeled on the basis of either Newtonian deterministic dynamics or a Langevin-type
stochastic dynamics, given their initial positions and velocities.

As the computer gets more powerful today and more accurate interatomic
potentials were developed, modeling of more realistic and more complicated systems
becomes possible. General materials always contain defects, such as grain boundary,
dislocation, cracks, void, vacancy, impurities, etc. To study these materials requires to
model ever-increasing system scale, at least millions of atoms. The massively parallel
(MP) computing hardware got improved in the last 10 years, so did the parallel

algorithms for MD simulations.

Chapter 2

2.2 Molecular dynamics (MD) simulations

Molecular dynamics (MD) is a kind of computer simulation techniques solving
the Newton's equations of motion with the time evolution for a collection of atoms
interacting via a potential U. The equations of motion for all atoms in a system are
integrated numerically by various finite differential methods at every time step. MD
simulations could generate detailed phase space information for a system, such as atomic
position and velocities at each time step. This information is also called the trajectory of
the system. Further analysis of the trajectory from MD simulations provides the linkage
between microscopic properties and average thermodynamic properties, such as pressure,
temperature, internal energy, etc.

There are several different forms of molecular dynamics to simulate different
ensembles. The original form of molecular dynamics generates the microcanonical
ensemble, or constant volume and constant total energy dynamics (NVE). Nose 2 added an
extra degree of freedom to describe the thermal bath behavior, such that the temperature
of the system will fluctuate with respect to the thermal bath temperature, this method can
achieve canonical ensemble, or constant volume and constant temperature dynamics
(NVT). Hoover3 further developed the Nose method to make the NVT calculation

simpler. Andersen developed a procedure to carry out isobaric-isoenthalpic ensemble, or
constant pressure and constant enthalpy dynamics (HPN), by making volume a dynamical
variable.
Parrinello and Rahman 5 generalized Andersen method to allow the changes in the
size and shape of the simulation cell. They define a new matrix h by h=(a,b,c), where a, b

Chapter 2

and c are the three vectors spanning the periodic repeating parallelepiped simulation cell.
In Parrinello and Rahman's theory, h becomes a dynamical variable to describe the shape
and size changes of simulation cell. The introduction of h into MD simulations makes it
possible to give a full description of the elastic properties of the system. Thus, one can
define the strain and stress tensor to be a new pair of extensive and intensive variables as
V and P for a thermodynamics system, which lead to constant thermodynamic tension
and constant enthalpy dynamics (HtN), or isobaric-isoenthalpic ensemble. The
introduction of h also clarified that the original NVE dynamics is actually NhE, and
original Nose constant NVT MD is actually NhT form of molecular dynamics, in which
the simulation box is kept unchanged not only in size but also in shape.
Ray and Rahman 6 have presented a detailed form of TtN dynamics, which
combines Nose constant-temperature theory with Parrinello-Rahman variable shape-size
form of molecular dynamics. One can also combine Nose's theory and Andersen's
changing volume dynamics to achieve constant TPN dynamics. The TPN dynamics is
suitable for isotropic liquid and gas phases, while the TtN dynamics can simulate elastic
deformation of solid state.

In the remainder of this section, we discuss MD simulation methods for four
ensembles (EhN, ThN, HtN and TtN). The Ray and Rahman's single Hamiltonian
formulation is used to cover all of these different forms of MD.

The Hamiltonian for the TtN form of MD has the form

H(s,1r,h,IT,f,P) =

7r;
Tr(IlIT)
L; if;G+u +---+ V Tr(t£) +-+ (3N + l)K T ln(f) ,(1)
2m;J
2W
2M

Chapter 2

where (si, lii) are the scaled coordinates and conjugate momentums of particle i, (h,TI) are
the coordinates and momentums of the simulation cell, and (f, P) are the Nose mass
scaling variable and its conjugate momentum. U is the potential energy, which is a
function of the position of atoms. The constants W ("piston mass") and M ("thermal
inertia") are parameters to make h and f satisfy dynamical equations. The tilde indicates
matrix transpose. T O is the thermal reservoir temperature, £ is the strain matrix which is
given by
(2)

where G is the metric tensor, G = hh, and ho is the reference state of the cell matrix h at
zero tension. VO is the reference volume, calculated from V0=det(h 0).
The usage of h matrix maps the simulation cell with any shape into a unit cell.
Thus the position and momentum (ri, pi) of physical particles are related to the scaled
particle variables (si, lii) by r; = hs; and P; = h -11l; If, and si range from O~ 1. Therefore,
the particle kinetic energy is represented by the first term in Hamiltonian (if we define the
physical momentum of the particle as P; = m; fhs;, then KE=

L P; I 2m; ), and the first

two terms in Eq. (1) are the Hamiltonian for N particles in the simulation system.
The elastic energy of the system given in the 4th term in Eq. (1), and the 3rd term
in Eq. (1) is similarly to the kinetic energy with the momentum of the h matrix. The last
two terms are a similar kinetic term and the potential term for f (the mass scale variable to
achieve constant temperature dynamics).

The equations of motion derived from Hamiltonian in Eq. (1) has the form of

Chapter 2

mi f 2··
Sia

(au !dr)s

!Ja - - -'l- -

r;,i

mi (j 2G-'G + 2fjPf')
J Sia,

(3a)

(3b)
Mj = _2K_E
_ _ _(3_N_+_I)_k8_T_

(3c)

where Pa~ is the microscopic stress tensor and the second term in Eq. (3b) is related to the
applied tension to the system. It is the difference of the system tension and applied
tension that causes the fluctuation of h matrix.

With no constraints, TtN dynamics requires the solutions of 3N+9+ 1 equation of
motions. (N is the number of movable particles with 3 degrees of translation freedom, h
matrix has 9 independent components, and one more degree of motion of j). We can get
the equations of motion for the other three dynamics from TtN dynamics by exerting
constraints. If the Nose variable f satisfies f = 0, f =I, then only the Eqs. (3a) and (3b)
are needed to be solved. This way, we reduce the constant TtN MD to the constant HtN
dynamics. Similarly, constraints on h as h = 0, h=constant lead to the equations of
motion reduced to the combination of Eqs. (3a) and (3c ), such that a constant ThN
dynamics is achieved. If f and h satisfy that j = 0, f = I and /2 = 0, h=constant, the only
equation of motion is Eq. (3a), which gives the constant EhN dynamics.

Chapter 2

2.3 Embedded atoms model (EAM) force fields (FF)

2.3.1 Physical foundation

The embedded-atom method (EAM) force field is a many-body potential for
computing the total energy of metallic systems, in which coordinate-dependent (or manybody) interactions are prominent7. In contrast, much simpler pair potentials always lead
to elastic constants C12=C44 (Cauchy relation) in cubic solids and the ratio of the vacancy
formation energy to cohesive energy as unity, which strongly deviate from the
fundamental properties of metallic solids 8 •
9 10

Daw and Baskes •

first proposed the EAM potential. They view the energy of

the metal as the energy obtained by embedding an atom into the local electron density
provided by the other atoms of the system. In addition, there is an electrostatic
interaction. The formula they used is

(4)

where G is the embedding energy defined as the interaction of the atom with the
background electron gas. The background electron density for each atom in the equation
is determined by evaluating at its nucleus the superposition of atomic-density tails from
other atoms. pa is the spherically averaged atomic electron density and U is an
electrostatic, two-atom interaction. A particular appealing aspect of the above EAM is its

Chapter 2

physical picture of metallic bonding, i.e., each atom is embedded in a host electron gas
created by its neighboring atoms.

Next, I will show how to derive the approximate expression as Eq. (4) for the
cohesive energy of a metallic system that is an explicit function of the positions of the
atoms.

The density functional expression for the cohesive energy of a solid 1s as
follows 11

coh

"f Z;P(!_)
dr +_!_ff PCiDp(rz) dr,dr - E

= G[p] + _!_",
Z;Zj L, R
L,;
2 ,,1

- _
1r

; 1

(5)

atom, '

r12

where the sums over i and j are over the nuclei of the solid, the primed sum indicates the
omission of the i=j term, Zi and i{ are the charge and position of the ith nucleus, the
integrals are over r (or ~ and r2 ), and r12=l ~ - r2 J. Eatoms is the collective energy of the
isolated atoms. G[p] is the kinetic, exchange, and correlation energy functional.

To go from Eq. (5) to Eq. (4), the following two assumptions are made.

(a) G[p] can be described by G[p] = g(p(r), Vp(r), V 2 p(r), ... )dr, where g is
the density and is assumed to be a function of the local electron density and its
lower derivatives,
(b) The electron density of the solid can be described as a linear superposition of
the densities of the individual atoms p_, (r) =

L Pt (r - i{).

Chapter 2

The first approximation is motivated by studies of the response function of the nearly
uniform electron gas. The second approximation is justified by the observation that, in
many metals, the electron distribution in the solid is closely represented by a
superposition of atomic densities. In addition, due to the variational nature of the energy
functional, errors in the assumed density should only affect the energy to second order. It
is also useful to define the embedding energy for an atom in an electron gas of some
constant density 75 (neutralized by a positive background):
G; [75;] = G[pt + 75;] - G[pt] - G[75;]. Using the above assumptions and the definition

for the embedding energy, the Eq. (6) can be obtained.

(6)

The error (Eerr) is a function of the background density 75;. Setting the error to
zero gives an equation for the optimal background density. The solution to Eerr=O is
discussed in detail by Daw 12 .

The EAM method has been applied successfully to study bulk and interface

problems, such as phonons 3, thermodynamics functions and melting point 14· 15 , liquid
metals

defects

17 19
- ,

grain boundary structure21 -25 , alloys 18 · 19 ·26 ·27 , segregation to grain

. 21-29 , mter
d"ff
. m
. all oys 30 '31 , an d fracture an d mec h amca
. l properties
. 32-38 .
b oun d anes
1 us10n
The EAM has been also applied to problems in surface structure 17 -I 9 ,39 -42 , adsorbate phase
43-47
· . 40 ,
d iagrams
, segregat10n
to surfaces 48-54 , surface structuraI ord er- d.1sord er trans1t10ns
surface ordered alloys

41 53
55 56
57 58
' , surface phonons ' , and clusters on surfaces · .

Chapter 2

Despite much success, EAM method will not work as well in the following two
cases: (1) where directional bonding is important, such as semiconductors and elements
from the middle of the transition series 59 and (2) where the Fermi-surface or bandstructure effects are important.

2.3.2 The qEAM FF

We develop and use the qEAM many-body Embedded-Atom-Model (EAM) type
force field (FF) for Ta. This FF is based on ab initio QM calculations and has been used
previously for molecular dynamics (MD) studies of the melting temperature of Ta as a
function of pressure6°, where it predicts values in excellent agreement with experiment. It
has also been used to characterize the nature of spall failure 61 .

The qEAM FF uses a functional form similar to that proposed by Chantasiriwan
and Milstein 62 . The total energy of system with atomic positions {rd is given by
E= LF(p;)+ L¢Cru),

(7)

with

P; = LfCru),

(8)

i#cj

where F(p) is the embedding energy, Pi is the total "electronic density" at site i, f( ru) is
the electron density function, between atoms i and j.
The electronic density is given by

Chapter 2

f(r)=

[l + a 1 cos(ar IV 3 ) + a 2 sin(ar IV 3 )

(9)

where Vis the volume per atom, a 1=0.07293238, a 2=0.15781672, a{llA)=21.79609053,
and /1=7.79329426.

The pair potential

where r,n(A)=4.81253968 is the cutoff radius. The parameters bi have the units of

eVIA

(4+i)

with ho= 6.50281587, b1 = -11.26455130, b 2 = 8.01451544, b3 = -2.97299223, b4

= 0.60004206, bs= -0.06222106, b6 = 0.00258801, and b 7 = -0.00000504.

The embedding function F(p) is determined from the Rose universal equation of
state 63 ,

1"
F(p) = EE0 5 (a * ) - L//J(r),

(lla)

* k a .3 +J4a *4) e -a' ,
EEOS(a *) =- E coh (l +a+

(1 lb)

where

with
(llc)

The parameters

entering the

definition

of the embedding energy are

aa(A)=3.32389219, Eco1z(eV)=8.154204, /4=0.207828, k=-0.00717801, and
f4= -0.00000504.

Chapter 2

2.3.3 Validation of the qEAM FF60

As input data to fit the qEAM FF, we use the following results from the linearized
augmented plane wave (LAPW) method with the generalized gradient approximation
(GGA) QM calculations.
(1) The zero temperature equation of state (EOS) of Ta for bee, fee, and hep
crystal structures for pressures up to - 500 GPa,
(2) the elastic constants,
(3) the volume relaxed vacancy formation energy also as a function of pressure,
(4) the equation of state for the A15 structure of Ta,
(5) the (100) surface energy in the bee Ta,
(6) the energies for shear twinning of the bee Ta.
To describe the properties of dislocation in bee Ta accurately, the qEAM FF must
reproduce the quantum calculation results for important quantities, such as equation of
state, elastic constants and energetics of homogeneously shear for bee Ta crystal. In the
following we show a detailed comparison between the qEAM FF and the data it was
fitted to.

A. Equation of state (EOS) for bee Ta

Figure 2-1 shows the energy [in Figure 2-l(a)] and pressure [in Figure 2-l(b)] as a
function of volume for bee Ta at T=0 K. The circles denote the LAPW GGA results and
the lines the qEAM FF calculations. The QM and the qEAM FF results agree to each
other. We also calculate the T=330 K EOS for bee Ta with the qEAM FF using

Chapter 2

isothermal-isobaric (NPT) MD with a Hoover3 thermostat and Rahman-Parrinello 64
barostat.

(a)

Ta
bccEOS

,.-..

-2

qEAM FF (lme)
QM-LAPW GGA (circles)

,>-,

--4

II)

---6

-8
-10

volume (A)

(b)

500
400

bccEOS
,.-..
ell

P-,

qEAM FF (line)
QM-1.APW GGA (circles)

300

200

0..
100

volume (A)

Chapter2

Figure 2-1. Zero temperature EOS for bee Ta, LAPW GGA and qEAM FF results: energy
[Figure 2-l(a)] and pressure [Figure 2-l(b)] as a function of volume. Circles denote
LAPW GGA results and lines show qEAM FF results.

In Table 2-1 we show the zero pressure volume, bulk modulus and its first

derivative with respect to volume for bee Ta at T=0K and 300K; we also show recent
compressibility data65 obtained in a diamond-anvil cell at room temperature. It is clear
that the qEAM FF reproduces the EOS for bee Ta very well.
Table 2-1. EOS parameters for bee Tantalum.

Vo (A 3)

BT (GPa)

Cn
(GPa)

C12
(GPa)

C44
(GPa)

LAPW-GGA

18.33

188.27

4.08

245.18

159.8

67.58

qEAMFF

18.36

183.04

4.16

272.54

137.57

69.63

FP LMTO GGA SC a

17.68

203

281

163

qEAMFF

18.4

176

4.9

Experiment b

18.04

194.7±4.8

3.4

264

159.7

82.2

T=0K

T=300 K

a Reference 66.
b Reference 65.

Chapter 2

B. Elastic constants

By fitting our energy-volume data [Figure 2-l(a)] to Rose's universal equation of
state, we obtained zero pressure volume (V 0), zero temperature bulk modulus (BT), and
its derivative with respect to pressure (B/) and, furthermore, the bulk modulus at
different pressures.

Static elastic constants [Cs = (C 11 -Cn)/2 and C 44] were obtained from strain
energies by straining the bee cell with volume conserving tetragonal and orthorhombic.
We calculate Cs using tetragonal strain of the cubic bee lattice:

ci = a(l + E,0,0),
b = a(O,l + E,0) ,

c = a(O,O,l/(l + c) 2 ) ,

(12)

where a is the cubic lattice constant of the system, a, b and care the lattice vectors and
E is the strain. Cs is related to the quadratic term of the strain energy

(13)

where E 0 is the energy of the unstrained system and VO is its volume.

Similarly C44 is obtained from the orthorhombic strain:

ci = a(l,E,0),
b = a(E,l,O),

c = a(O,O,l/(l - £

)) ,

(14)

the shear constant C44 is obtained from
(15)

Chapter 2

Table 2-1 reported the obtained bulk modulus and elastic constants (C 11 , C 12 and
C44 ) at 0 K and 300 K for bee Ta using the qEAM FF and QM calculations. Figure 2-2
shows the elastic constants [bulk modulus BT=(C 11 +2C 12 )/3, Cs=(Cu-C 12 )/2 and C44 ] as a
function of pressure obtained with the qEAM FF (filled circles and full lines) and the
LAPW results. While the agreement in BT is excellent and that for Cs is good, the qEAM
FF greatly underestimates C44 for high pressures.

Ta
elastic constants
1500

QM (LAPW) dashed lines & open symbols
qEAM full lines & symbols

Bulk Modulus (circles)
c 44 (diamonds)
500

c 5 (squares)

400

600

pressure (GPa)

Figure 2-2. Zero temperature elastic constants for Ta, LAPW GGA and qEAM FF results.
Circles show bulk modulus [(C 11 + 2C 12 )/3]; diamonds show C44 and squares represent Cs
= (C11-C!2)/2. qEAM FF results are shown with filled symbols and full lines and ab initio
LAPW results with open symbols and dashed lines.

Chapter 2

C. Energetics of homogeneously sheared bee crystal

The ideal shear strength is defined to be the stress separating elastic and plastic
deformation when a homogeneous shear is applied to a perfect crystal. It gives an upper
bound for the shear strength of the material. The shear transformation is in the direction
of the observed twining mode and deforms the crystal into itself66 • 67 . For bee crystal we
use the following transformation of the cell vectors 66 • 67

1 -

a=-[111]+ r;;:;-[111],
-vl8

- 1 s -b = -[111] + r;;:;- [l 11],
-vl8
(16)

when the shear variable sis equal to the twinning shears = Stw = 2 -1/Z the lattice vectors
[a= 1/3 [ 212], b=l/3 [122] and c = [11 1 ]] form a bee structure, twin of the initial one.

In this way one can calculate the energy along the shear path,
W(s) =e(V, s) - e(V, s =0),

(17)

where e(V, s) is the energy per atom of the deformed system and e(V, s = 0) is the perfect
crystal energy. The energy barrier associated with this transformation is W max = W (s =
0.5). The corresponding stress is defined as:
i-(s) = 1 dW(s),
V ds

(18)

The ideal shear strength (Tmax) is defined as the maximum stress along the path. Figure 23 shows energy and stress as a function of shear using the qEAM FF for zero pressure
volume.

Chapter 2

(a) o. 25 ,---y----r--r---.--r--.--.----.---~--.-.----.--.-,--,r-r--.--r--.---.--.--T--,
0.2

> 0.15

-t
l1J

@ 0.1

qEA1v.1
V = 18.36A

w max

6 0.05
-0. 05 L......J.......I......L..-..__._.....1..._.__..__.__...__.__....__.,__.__.__....___.__.__,__....___.__._..............

--0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.1

shear x/stw

(b)

10
Ta

qEA1v.1
V = 18.36 A

ti.I
ti.I

tal1J

-5

-10'--'---'----'--'-..........---'----'--'-...........---'---'---'-..........___._--'---'-..........---'-...............

--0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.1

shear x/stw

Figure 2-3. Ideal shear strength of Ta using qEAM FF at zero temperature and volume

V=18.36 A • We show energy W(s) [Figure 2-3(a)] and stress i-(s) [Figure 2-3(b)] as a
function of shear.

Soderlind and Moriarty calculated W(s) and i-(s) for Ta at different volumes, from
first principles. In developing the qEAM FF we used Wmax for V=17.6186 A3 and

V=I0.909 A as part of the training set. In Table 2-2 we show a comparison between the

Chapter 2

first principles results and the ones obtained using the qEAM FF. We can see that the
qEAM results 66 are in very good agreement with the ab initio calculations.
Table 2-2. Shear deformation in the observed twinning mode in Ta.

Volume (A 3)

Wmax (eV)

'tmax (GPa)

Wmax (eV)

'tmax (GPa)

FP LMTO GGA SC a

qEAMFF
18.360000

0.188

7.14

17.618602

0.2

8.0

0.194

7.37

15.143996

0.26

12.05

0.276

12.4

10.9090116

0.43

28.2

0.566

36.2

a Reference 66.

Chapter 2

2.4 References

1. M. P. Allen and D. T. Tildesley, Computer Simulation of Liquids, (Oxford, 1987).
2. S. Nose, J. Chem. Phys., 81, 511 (1984).
3. W. G. Hoover, Phys. Rev. A, 31, 1695 (1985).
4. H. C. Andersen, J. Chem. Phys., 72, 2384 (1980).
5. M. Parrinello and A. Rahman, Phys. Rev. Lett., 45, 1196 (1980).
6. J. R. Ray and A. Rahman, J. Chem. Phys., 82, 4243 (1985).
7. M. S. Daw, S. M. Foiles and M. I. Baskes, Mater. Sci. Rep., 9,251 (1993).
8. F. Ercolessi, M. Parrinello and E. Tosatti, Philos. Mag. A, 58, 213 (1988).
9. M. S. Daw and M. I. Baskes, Phys. Rev. Lett., 50, 1285 (1983).
10. M. S. Daw and M. I. Baskes, Phys. Rev. B, 29, 6443 (1984).
11. P. Hohenberg and W. Kohn, Phys. Rev. B, 37, 10411 (1964).
12. M. S. Daw, Phys. Rev. B, 39, 7441 (1989).
13. M. S. Daw and R. D. Hatcher, Solid State Commun., 56,697 (1985).
14. S. M. Foiles and M. S. Daw, Phys. Rev. B, 38, 12643 (1988).
15. S. M. Foiles and J.B. Admas, Phys. Rev. B, 40, 5909 (1989).
16. S. M. Foiles, Phys. Rev. B, 32, 3409 (1985).
17. M. S. Daw and M. I. Baskes, Phys. Rev. B, 29, 6443 (1984).
18. S. M. Foiles, M. I. Baskes and M. S. Daw, Phys. Rev. B, 33, 7983 (1986).
19. S. M. Foiles, M. I. Baskes and M. S. Daw, Phys. Rev. B, 37, 10387 (1988).

Chapter 2

20. S. M. Foiles, M. I. Baskes and M. S. Daw, in Interfacial Structure, Properties and

Design, eds. M. H. Yoo, W. A. T. Clark and C. L. Briant, (Materials Research
Society, Pittsburgh, PA, 1988).
21. S. M. Foiles, Acta Metall., 37, 2815 (1989).
22. M. J. Mills, G. J. Thomas, M. S. Daw and F. Cosandey, Mater. Res. Soc. Symp.
Proc., 159, 365 (1989).
23. U. Dahmen, C. J. D. Hetherington, M. A. O'Keefe, K. H. Westmaccott, M. J. Mills,
M. S. Daw and V. Vitek, Philos. Mag. Lett., 62, 327 (1990).
24. M. J. Mills and M. S. Daw, Mater. Res. Soc. Proc., 183, 15 (1990).
25. M. J. Mills, M. S. Daw, G. J. Thomas and F. Cosandey, Ultramicroscopy, 40, 247
(1992).
26. S. M. Foiles and M. S. Daw, J. Met., 39, 39 (1987).
27. S. M. Foiles, in High-Temperature Ordered Intermetallic Alloys, Vol. 2, eds. N.S.
Stoloff, C. C. Koch, C. T. Liu and 0. Izumi (Materials Research Society, Pittsburgh,
PA, 1987).
28. S. M. Foiles, Phys. Rev. B, 40, 11502 (1989).
29. A. Seki, D. N. Seidman, Y. Oh and S. M. Foiles, Acta Metall., 39, 3167 (1991).
30. J.B. Admas, S. M. Foiles and W. G. Wolfer, J. Mater. Res., 4, 102 (1989).
31. J. B. Admas, S. M. Foiles and W. G. Wolfer, in Atomistic Simulation of Materials,
eds. D. Srolovitz and V. Vitek, (Plenum Press, New York, 1989).
32. M. I. Baskes, S. M. Foiles and M. S. Daw, J. Phys., CS, 483 (1988).

Chapter 2

33. M. S. Daw, M. I. Baskes, C. L. Bisson and W. G. Wolfer, in Modeling Environmental

Effects on Crack Growth Processes, Metallurgical Society Fall Meeting, Toronto,
Canada (1985), eds. R.H. Jones and W.W. Gerberich.
34. M. S. Daw and M. I. Baskes, in NATO Advanced Workshop on Physics and

Chemistry of Fracture, Bad Reichenhall, Weat Germany, eds. R. M. Latanision and
R.H. Jones (Nijhoff, Dordrecht, 1986).
35. M. I. Baskes and M. S. Daw, in Fourth International Conference on the Effect of

Hydrogen on the Behavior of Materials, Jackson Lake Lodge, Moran, WY, eds. N.
Moody and A. Thompson (The Minerals, Metals, and Materials Society, Warrendale,
PA, 1989).
36. R. G. Hoagland, M. I. Baskes, M. S. Daw and S. M. Foiles, J. Mater. Res., 5, 313
(1990).
37. R. G. Hoagland, M. S. Daw, S. M. Foiles and M. I. Baskes, in Atomic Scale

Calculations of Structure in Materials, eds. M. S. Daw and M.A. Schluter (Materials
Research Society, Pittsburgh, PA, 1990).
38. R. G. Hoagland, M. S. Daw and J. P Hirth, J. Mater. Res., 6, 2565 (1991).
39. S. M. Foiles, Surface Sci., 191, L779 (1987).
40. M. S. Daw and S. M. Foiles, Phys. Rev. Lett., 59, 2756 (1987).
41. M. S. Daw and S. M. Foiles, in 2 nd International Conference on the Structure of

Surfaces, Amsterdam, The Netherlands, eds. J. F. van der Veen and M.A. Van Hove,
(Springer, Berlin, 1987).
42. M. S. Daw and S. M. Foiles, J. Vacuum Sci. Technol. A, 4, 1412 (1986).
43. S. M. Foiles and M. S. Daw, J. Vacuum Sci. Technol. A, 3, 1565 (1985).

Chapter 2

44. M. S. Daw and S. M. Foiles, in 1st International Conference on the Structure of

Surfaces, Berkeley, CA, eds. M. A. Van Hove and S. Y. Tong (Springer, Berlin,
1985).
45. T. E. Felter, S. M. Foiles, M. S. Daw and R. H. Stulen, Surface Sci., 171, L379
(1986).
46. M. S. Daw and S. M. Foiles, Phys. Rev. B, 35, 2128 (1987).
47. T. L. Einstein, M. S. Daw and S. M. Foiles, Surface Sci., 227, 114 (1990).
48. S. M. Foiles, Phys. Rev. B, 32, 7685 (1985).
49. S. M. Foiles, J. Vacuum Sci. Technol. A, 4, 761 (1986).
50. S. M. Foiles, in Computer-Based Microscopic Description of the Structure and

Properties of Materials, eds. J. Broughton, W. Krakow and S. T. Pantelides
(Materials Research Society, Pittsburgh, PA, 1986).
51. S. M. Foiles, in Physical and Chemical Properties of Thin Metal Overlayers and
Alloy Surfaces, eds. D. M. Zehner and D. W. Goodman, (Materials Research Society,
Pittsburgh, PA, 1987).
52. S. M. Foiles, J. Vacuum Sci. Technol. A, 5, 889 (1987).
53. S. M. Foiles, Surface Sci., 191, 329 (1987).
54. S. M. Foiles, in Surface Segregation and Related Phenomena, eds. P.A. Dowben and
A. Miller (CRC Press, Boca Rotan, 1990).
55. J. S. Nelson, M. S. Daw and E. C. Sowa, Phys. Rev. B, 40, 1465 (1989).
56. J. S. Nelson, E. C. Sowa and M. S. Daw, Phys. Rev. Lett., 61, 1977 (1988).
57. P. Schwoebel, S. M. Foiles and G. Kellogg, Phys. Rev. B, 40, 10639 (1989).
58. A. F. Wright, M. S. Daw and C. Y. Fong, Phys. Rev. B, 42, 9409 (1990).

Chapter 2

59. A. E. Carlsson, in Solid State Physics: Advances in Research and Applications, eds.
H. Ehrenreich and D. Turnbull, (Academic Press, Boston, 1990), Vol 43, Pl.
60. A. Strachan, T. <:;agm, 0. Gtilseren, S. Mukherjee, R. E. Cohen, and W. A. Goddard,

First Principles Force Field for Metallic Tantalum, Phys. Rev. B, submitted.
61. A. Strachan, T. <:;agm and W. A. Goddard, Phys. Rev. B, 63, 060103 (2001).
62. S. Chantasiriwan and F. Milstein, Phys. Rev. B, 53, 14080 (1996).
63. J. H. Rose, J. Ferrante, and J. R. Smith, Phys. Rev. Lett., 47, 675 (1981).
64. M. Parinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
65. H. Cynn and C. Yoo, Phys. Rev. B, 59, 8526 (1999).
66. P. Soderlind and J. A. Moriarty, Phys. Rev. B, 57, 10340 (1998).
67. A. T. Paxton, P. Gumbsch and M. Methfessel, Philos. Mag. Lett., 63, 267 (1991).

Chapter]

Chapter 3 Core structure and core energy of 1/2a screw dislocation
3 .1 Overview
The core energy of an isolated dislocation is an essential parameter in modem
plasticity theory1. Xu et al. 2 , Beigi and Arias 3 , and Yang et al. 4 have computed the core
energy of the a/2<111> screw dislocation for Mo and Ta by applying anisotropic elastic
theory. There are two important parameters in their calculations.
1. Anisotropic shear modulus. It can be derived directly from the elastic constants of

the perfect crystal 5 or extracted from the atomistic simulations of the dislocation
energy as a function of cell size (as will be shown below). A good agreement between
the results from two approaches indicates a quantitative correspondence between the
elasticity theory and atomistic simulations.
2. Dislocation core radius. In previous studies2·3 , an approximate value 2b was obtained
by fitting the strain energy for cylinders containing a dislocation with various radii to
anisotropic elastic theory. Since the computed core energy depends strongly on the
core radius, a physically based definition of the dislocation core radius is necessary.
To provide an increased atomistic insight into the nature of the dislocation core,
we computed the dislocation core energy using two approaches (an atomistic model and a
continuum model). By comparing and contrasting these two approaches, we obtained a
consistent core energy and core radius of the a/2<111> screw dislocation in Ta.

Chapter3

3.2 Dislocation core structure
3 .2.1 Construction of the dislocation quadrupole
To investigate the core structure and determine the core energy of l/2a
screw dislocations in Ta, we use a quadrupole arrangement of dislocations in a periodic
simulation cell. Thus, two of the dislocations have Burgers vector h=a/2[111] and the
other two have h=a/2[-1-1-1]. This arrangement (Figure 3-1) leads to little positional
misfit of atoms across the cell boundary due to the effect of periodic images.

[-110]

3/4

[111]

[11-2]

1/4

3/4

Figure 3-1. The geometrical arrangement of the dislocation quadrupole simulation cell.
The [111], [-110] and [11-2] boundaries of simulation cells are periodic boundaries.

The initial configuration of the dislocation quadrupole was constructed using
isotropic elastic theory. The displacement of each atom along the [111] direction is given
by

(1)

Chapter3

where the sum runs over all dislocations on the

plane.

(f)i

1s the

counterclockwise angle on the (111) plane from the [11-2] direction to the vector joining
the center of dislocation I to the atomj, while b1 is the Burgers vector of dislocation I and
b1 is equal to la/2<111>1=2.88 A.

dislocations are initially located at the geometric centers of a triangle
surrounded with three [111] columns of atoms. For the bee structure, there are two kinds
of dislocation core configurations that can be transformed to each other by reversing the
Burgers vector; they are called "easy core" [Figure 3-2(c)] and "hard core" [Figure 32(b)]2. The easy core is the low energy form and the only one we find from
minimizations. Indeed, the dynamical simulations (Chapter 4) show that the dislocation
moves from one easy form to an adjacent one avoiding the high-energy hard core.

A BCABC

ABCBCA

ABCCAB

(a)

(b)

(c)

3-2. The [-110] view
sequence of {1

} planes on

the a/2<111> screw dislocation
sides of

Ta. The stacking

dislocation is shown as (a) 'ABCABC'

Chapter]

bulk bee crystal, (b) 'ABCBCA' in the "hard core" screw dislocation and (c) 'ABCCAB'
in the "easy core" screw dislocation.

3.2.2 The differential displacement (DD) map

[111)
(11-2)

•> • > o~ •

t,,

• ..-:=o -..: • <• .co

(1-1 OJ

O••O

"a•

•O••O

■ < •

•..::;o<•..:::•..:::o

•O••O

•O••O••O••

•~">a>•

••a

Figure 3-3. Differential displacement (DD) map for the equilibrated dislocation
quadrupole in which there are 5670 atoms and each dislocation is 7 b long.

Some important features of dislocation cores can be visualized using differential
displacement (DD) maps 6 . Figure 3-3 shows the DD map for a quadrupolar arrangement
of dislocations after relaxing the atomic positions (using the qEAM FF). In this map, the
atoms are represented by circles and projected on a (111) plane of the bee lattice. The

Chapter3

arrows in the DD map indicate the relative displacements of neighboring atoms in the
[111] direction with respect to their positions in the perfect bee crystal. The direction of
the arrow represents the sign of the displacement and the magnitude of the arrow is
proportional to the relative displacement between the corresponding atoms. When an
arrow spans the full distance between two atoms, the relative displacement is b/3.

Figure 3-3 shows that the equilibrium dislocation cores have threefold symmetry
and spread out in three <112> directions on {110} planes in the DD map. There are 6
equivalent <112> directions on the (111) plane, so there are two kinds of "easy core"
configurations, each with the threefold symmetry. The dislocations with different core
configurations are energetically degenerate both in terms of core energy and in terms of
elastic energy. In a dislocation quadrupole, changing the dislocation core configuration
does not affect the equilibrium energy of the simulation cell.

3.2.3 Polarization of the dislocation

In addition to the relative displacement of neighboring atoms in the [111]
direction, we computed the difference between the atomistic relaxation and the
predictions of isotropic elasticity theory [Eq. (1)] for the atomic displacements parallel to
the Burgers vector. We find that except for the 6 columns of atoms closest to the
dislocation line, elastic theory and the atomistic relaxation lead to atomic displacement
differences within (-0.05 A, 0.05 A). This shows that the continuum theory accurately
describes the elastic displacement field of a screw dislocation, failing only for the
innermost 6 columns of atoms within the dislocation core.

Chapter3

An interesting feature from the atomistic simulations is that after relaxation the
three central columns of atoms of the dislocation translate simultaneously by 0.267 A
(-0.09 b) either in the [111] direction or the [-1-1-1] direction. This phenomenon is called

dislocation polarization7. Dislocation polarization can be quantified by the simultaneous
displacement in the [111] direction of the central three atoms at the dislocation core4 . The
[111] polarization leads to the dislocation core spreading along the [-1 -1 2], [-1 2 -1],
and [2 -1 -1] directions in DD map. For the opposite polarization (in the [-1-1-1]
direction) the dislocation cores spread out along the [1 1 -2], [1 -2 1], and [-2 -1 -1]
directions in DD map.

The two dislocations on the left of Figure 3-3 have [-1-1-1]

polarization while the two dislocations on the right have [111] polarization. The above
relationship between the dislocation polarization and the directions in which the
dislocation core spreads is independent of the orientation of the Burgers vector.

Figure 3-4 shows differential displacement maps and relaxation maps for the four
types of 1/2a screw dislocations: N+ (polarized in the [-1-1-1] direction and with
b=l/2a[ll 1]) dislocation, P+ (polarized in the [111] direction and with b=l/2a[l 11])
dislocation; N- (polarized in the [-1-1-1] direction and with b=l/2a[-1-l-1]) dislocation,
and P- (polarized in the [111] direction and with b=l/2a[-1-l-1]) dislocation.

Chapter3

i •
•o.!t!.o.

-0.267

0.123
-0.267

0.123

(b) P+

o •-•--o
I\

0 ~.m •

•,/~o •

l•

~.m 0

TO •

.. -· . • • r • •

0.2671

0.267
0.267
-0.123

()

• •

0.123
-0.267
-0.267

• 0 •

r-0.267

a'·"'•

(d) P-

• •

r•

0 •

j •'·"' o

Chapter3

Figure 3-4. The equilibrated dislocation core configurations for the 1/2a screw
dislocation in Ta. The circles represent the projected atoms in the (111) plane. The open,
shaded or black circles indicate that the atoms are in three consecutive (111) layers of bee
lattice. However, the arrows in two columns of figures have different meanings. The left
column shows the differential displacement map, in which the arrow indicates the
displacement in [ 111] direction (perpendicular to the map) of the neighboring atoms
relative to their positions in the perfect bee crystal. The direction of the arrow represents

the sign of the displacement and the magnitude is proportional to the relative
displacement between corresponding atoms. When the arrow touches the two atoms, the
relative displacement between these two atoms is 1/3 b. For clarity, the relative
displacements less than 1/12 b are not shown in the figure. The right column shows the
relaxation map, in which the arrow from each atom indicates the relaxation (parallel to
the dislocation line) relative to the displacement field predicted by isotropic elastic
theory. The magnitudes of such relaxation (in angstrom) for the central six columns of

atoms (the relaxation for the other atoms is less than 0.05A) are printed next to the
corresponding atom. Four types of energy degenerate dislocation core configurations are
distinguished in terms of the relaxation direction of the three central columns of atoms
(downward denoted as "N" and upward denoted as "P") and Burgers vector (a/2[111]
denoted as "+" and -a/2[111] denoted as "-").

Chapter3

3.2.4 Comparison to other calculations

Using the qEAM FF, we find a polarization of the dislocation core in which the
relative displacement field has the threefold symmetry, rather than the sixfold symmetry
of (111) planes of the bee lattice. However, we must be cautious about these results.
Previous studies for bee metals have led to both asymmetric (threefold symmetric) and
symmetric (sixfold symmetric) core structures for the a/2<111> screw dislocation.

Xu and Moriarty et al. 2 ·8·9 obtained an asymmetric core structure for bee Mo using
the multi-ion interatomic potentials from the model generalized pseudopotential
theory (MGPT).

Vitek arrived at an asymmetric core structure for bee crystal employing pair-wise
interatomic potentials from Ref. 6.

Experimentally, Sigle investigated the dissociation of the a/2<111> screw
dislocation in Mo using high-resolution transmission electron microscopy
(HRTEM), obtaining a broken symmetry consistent with an asymmetric
dislocation core configuration 10 •

Duesbery and Vitek 11 found a symmetric core for group VB metals (V, Nb, and
Ta) but an asymmetric core for group VIB metals (Cr, Mo, and W) using FinnisSinclair (F-S) type central-force many body potentials.

Ab initio density functional theory calculations for both Mo and Ta led to a
symmetric dislocation core both with periodic boundary conditions 3 and lattice

Green's function boundary conditions 12 .

Chapter3

Yang et al. 4 showed the core structure of the screw dislocation for Ta is very
sensitive to the mechanical conditions and varies dramatically from symmetric to
asymmetric with increasing pressure.

Summarizing, the weight of evidence is in favor of a symmetric dislocation core, but the
evidence remains inconclusive.

The main difference between asymmetric [DD map in Figure 3-5(b)] and
symmetric [DD map in Figure 3-5(a)] core is the polarization. The symmetric core does
not have polarization. The extent of polarization for an asymmetric core is [-b/6,0) and
(0, b/6].

.... O-+-e

'• 0

(a)

• '~.'

• •

-o
I'

-- • - •
• • 0 • •

(b)

Figure 3-5. The DD maps for the equilibrium dislocation core structures. (a) The
dislocation polarization p is equal to 0.00 b and the corresponding core is symmetric,
while (b) the dislocation polarization p is equal to 0.09 b and the corresponding
dislocation core is asymmetric, spreading along three <112> directions on the (110)
planes.

Chapter3

Starting from the asymmetric core structure, we constructed a symmetric
dislocation core by translating the central six columns of atoms along the [111] direction
and then relaxing the full unit cell while fixing the positions of central six columns of
atoms in the [111] direction. The resultant symmetric core has atomic displacements
similar to those in the literature 3 · 11 ' 12 . We calculate using the qEAM FF that the
dislocation core energy per Burgers vector b for symmetric core is 1.440 eV, which is
0.040 eV higher (2.9 %) than the relaxed asymmetric core.

3.3 Dislocation core energy

3.3.1 Atomistic approach

We define the strain energy associated with each atom as in Eq. (2):

E; = F(p;) +

~ L. ¢(r;) - Ecoh'

(2)

j'#l

where g·oh is the atomic cohesive energy in perfect crystal.

The atomic strain energies calculated using Eq. (2) for the relaxed dislocation
quadrupole with 5670 atoms (with cell size: X=9<112>a, Y=15<110>a, and
Z=7/2<111>a) are displayed in Figure 3-6(a). Here each atom is projected on the (111)
plane and drawn as a circle whose radius is proportional to its atomic strain energy. Most
atoms have very small strain energy; only 12 atoms close to the dislocation line have
significant strain energies. Figure 3-6(b) shows the atomic strain energy distribution per
dislocation per b for the same dislocation quadrupole cell. Here we see that the atomic
strain energy of the six atoms close to the dislocation line is 0.15 eV to 0.17 eV, while

Chapter3

another 6 atoms have atomic strain energies ranging from 0.06 eV to 0.08 eV. These 12
atoms near the center of the dislocation in Figure 3-6(a) are denoted as A, B, C and D in
Figure 3-6(b), in decreasing atomic strain energy order. Except for these 12 atoms, all
other atoms have atomic strain energies less than 0.05 eV. Based on these observations,
we define the core of the dislocation to be formed by the 12 atoms with higher strain
energy per Burgers vector. This leads to the dislocation core energy of Ec=l.400 eV per
Burgers vector b .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .N.:...5670 .

A······ ····.,,B ···· ·

. . . . . . . . . . . . . . . . . . . . .

. . . . .

···-~BO·········· ·O·•· · · · ·
· · · o
····BAB•··········

AO·····
B · · · · • · ·

· · · · • · o · · · · · · · · · · · o

...................•......

. . . . . . . . . . . . . . . . . . . . . . .
...................•......

· · · ·•·O· · · · · · · · · · -~B- · · · · · ·
···-~AB•··········
O· · · · ·

· · · ·OA

· ·········•BA•·····

· · BO · · · · · · · · · · · o · • · · · · ·

.....•.....................

(a)

Chapter]

!Quadrupole, N=5670 I
elastic region

core region ( 12 atoms I

....,.,
C.

170.S

0.GJ 0.02 0.03 0.04 0.05 0.06 0.07 0.08 OJN 0.1 0.1 I 0.12 O.l.30.14 0.15 0.16 0.17

strain energy (eV)

(b)

Figure 3-6. (a). The <111> projection of atomic strain energy distribution for a lb thick
slab in an equilibrated dislocation quadrupole in which there are 5670 atoms and each
dislocation is 7b long. Atoms are represented as circles and the radius of circle is
proportional to the strain energy of the atom. (b ). Histogram of atomic strain energy
distribution for a 1b segment of the dislocation obtained from the same quadrupole
simulation. The number of atoms in each energy bin is shown on the top of the
corresponding bar.

In Figures 3-6 (a) and (b), the atoms with atomic strain energy ranging from 0.165
e V to 0.170 e V are denoted as 'A' and atoms denoted as 'B' have atomic strain energy
ranging from 0.156 eV to 0.157 eV. Atoms 'C' in Figure 3-6(b) are represented as black

Chapter3

circles in Figure 3-6(a), while atoms 'D' in Figure 3-6(b) are represented as white circles
in Figure 3-6(a). The atoms labeled by 'C' have higher strain energy than the atoms 'D'.

~- --0\

- ----

-----0

0---- -

(]

t3-----0' 0

Figure 3-7. Schematic drawings for a dislocation core in the atomistic model and in the
continuum model. Atoms are projected in (111) plane and represented by circles. The
black circles indicate the atoms constituting the dislocation core according to the
atomistic model, the other atoms are in the elastic region and drawn as white circles. The
dotted line connects the non-core atoms that most closely encircle the dislocation core
providing a cutoff boundary for the atomistic model. The solid circle whose radius is the
average distance from the dislocation center to the atomistic cutoff boundary (2.287 b) is
the dislocation core radius for the continuum model.

Chapter3

Corresponding to this definition, we define the core radius as the average distance
from the dislocation center to the closest non-core atoms encircling the core region. This
is shown in Figure 3-7 as the dotted line. This leads to a core radius of rc=2.287 b. Note
in Figure 3-7 that the 12 atoms of the core are not the 12 closest atoms to the dislocation
center.

Summarizing, the atomic strain energy distribution provides a criterion for
distinguishing which atoms are inside and outside the core region of the dislocation. This
leads to a core energy of Ec=l.400 eV/b and a core radius of rc=2.287 b for the l/2<111>
screw dislocation in Ta.

3.3.2 Continuum approach

The presence of dislocations leads to strain in the crystal. The total strain energy
can be considered as the summation of the dislocation core energy and the elastic energy
outside the dislocation core. The latter, including the dislocation self-energy and
interaction energy, can be calculated using elasticity theory. Inside the dislocation cores
the strains are too large for elasticity theory to apply.

The total strain energy per Burgers vector for two parallel straight dislocations
with equal and opposite Burgers vectors at a separation dis
E = 2Ec (rJ + 2Kb ln(-),

(3)

where r, is the core radius of the dislocation and the elastic modulus K can be expressed
as 5

Chapter3

K = {S11 l[S44(S11S44 -S1 s)]} 2

4n-

(4)

Sn, S44 , and S 15 are the modified elastic compliance constants. The details of

determining them from standard elastic constants of the cubic crystal can be found in
Ref.5.

Summing the pair interactions in Eq. (3) leads to the total energy per dislocation
per bin a dislocation quadrupole cell as in Eq. (5) 3

(5)

where d 1 and d2 are the distances between the dislocations along the [-11 O] and [11-2]
directions, A(d1/d2 ) is a convergent summation of all pair interaction effects and is related
to the geometry of the simulation cell. In Eq. (5), the core energy Ec(rc) and effective
elastic parameter Kb 3 are constants, leading to a total strain energy that varies linearly
with the scaled elastic energy [ln(d1lrc)+A(d1ld2)], as the size of simulation cell is
changed. Plotting the total strain energy versus the scaled elastic energy, we determine
the effective elastic modulus K from the slope and the core energy Ec(rc) from the
intercept. Obviously, the core energy Ec(rc) obtained in this way depends on the choice of
re, while K does not.

To determine Ec(rc) and K using Eq. (5), we simulated quadrupole arrays of
dislocations for various system sizes ranging from 1890 atoms (40.71A by 42.31A by
20.15A) to 51,030 atoms (219.8A by 211.5A by 20.15A) and optimized the atomic

Chapter3

coordinates by minimizing the energy. The geometrical parameters of simulation cells,
numbers of the atoms and the obtained strain energies are shown in Table 3-1.
Table 3-1. Table of size of simulation cells, number of atoms per simulation cell and
strain energy per dislocation per Burgers vector. X, Y and Z are the cell parameters for
the simulation cells. Xis in the unit of a[ll-2], Y is in the unit of a[-110] and Z is in the
unit of a/2(111]. IXI, IYI and IZI are the size of simulation cells in unit of A.

IXI (A)

IYI (A)

IZI (A)

N (/cell)

E (eV/b)

40.71

42.31

20.15

1890

1.833

40.71

51.71

20.15

2310

1.891

56.99

42.31

20.15

2646

1.927

56.99

51.71

20.15

3234

2.039

73.28

51.71

20.15

4158

2.095

73.28

70.51

20.15

5670

2.265

171.0

155.1

20.15

29106

2.912

219.8

211.5

20.15

51030

3.150

In Figure 3-8, we show that the total strain energy for various simulation cells as a
function of the scaled elastic energy [Zn( d1lrc)+A( d1/d2 )]. Our results show the linear
dependence expected from Eq. (5). Taking the core radius as 2.287 b obtained from the
atomistic method and using the linear fit of our data in Figure 3-8, we determine the
"easy core" dislocation to have the core energy of 1.404 eV/b.

This is in excellent

agreement with the dislocation core energy of 1.400 e Vlb calculated directly using the
atomistic model. This linear fit leads to an elastic modulus of K = 3.3497 x 10-2 eV/A 3 .
Alternatively using the computed elastic constants for the bee crystal from the qEAM FF,

Chapter3

we obtain the K = 3.3492 x 10-2 eV/A 3 , which is within 0.02% of the value derived from
fitting Eq. (5) in Figure 3-8.

3.5

-->
.c

2.5

C)
lr-

,c
«f

Ir-

+"'
(JJ

1.5

0.5

1 .5

2.5

scaled elastic energy (eV/Kb )
Figure 3-8. The strain energy (E/b) as a function of the scaled elastic energy
(ln[ d1lrcl +A(d1ld2)) obtained from the simulations of different sized dislocation

quadrupoles. The dotted line represents the linear fitting
Elb=l .4041 +O. 79899[ln( d1lrc)+A( d1ld2)] with rc=2.287 b leading to a dislocation core

energy of 1.4041 eV/b. The number of atoms in each simulation cell is specified in the
figure.
Beigi and Arias 3 calculated the core energy of l/2a <111> screw dislocation in Ta
from ab initio methods. They used r, = 2b which led to Ee= 0.86 eV/b. Using our force
fields, with the same core radius, we compute the core energy to be Ee = 1.297 eV/b.

Chapter3

Thus, the dislocation core energy based on ab initio QM 1s 34% lower than our
calculation.
In the MGPT calculations on Ta 4, a value of re= 1.75 b was used, leading to Ee=
0.60 eV/b. This can be compared to the value of Ec=l.190 eV/b that our results in Figure
3-8 would lead to for re= 1.75 b. Thus, the dislocation core energy based on the MGPT
force field is 50% lower than our calculation.

Table 3-2 compares the predicted lattice parameter and elastic constants for bee
Ta from the ab initio calculations 3,4· 12 · 14 , MGPT 4, and qEAM computations as well as the
experiments 13 . Predicted elastic constants differ by tens of a percent, especially for C44 . It
is likely that the predicted dislocation core energies from different computations would
have at least the same magnitude of diversity. However, even with such energy
differences, we expect that the character and properties of dislocations can be understood
properly from the simulations.

Table 3-2. Experimental and theoretical values of lattice parameter a (A), elastic
constants C11 (0Pa), C12 (0Pa) and C44 (0Pa). Errors with respect to the experimental
values for different ab initio pseudopotentials, MGPT and qEAM force field calculations
are shown.

Chapter]

Exp. a

a (A)

Cn (0Pa)

C12 (0Pa)

C44 (0Pa)

3.30

266

158

87.4

Ab Initio b

Value

3.25

304

182

(Beigi et al.)

Error

-1.5 %

14 %

15 %

-25 %

Ab lnitio C

Value

3.30

265

155

91.3

(Soderlind et al.)

Error

-0.4 %

-1.9 %

4.5 %

Ab lnitio d

Value

3.23

291

175

52.9

(Woodward et al.)

Error

-2.1 %

9.4 %

11%

-39 %

Ab Initio e

Value

3.36

244

160

66.3

(Gtilseren et al.)

Error

1.8 %

-8 %

1.3 %

-24 %

MGPTFFf

Value

3.30

266

161

82.5

(Yang et al.)

Error

1.9 %

-5.6 %

qEAMFF

Value

3.32

273

138

69.6

(present work)
a Reference 13.

Error

0.6%

2.6 %

-13 %

-20 %

Reference 3.

The

total-energy plane-wave density-functional pseudopotential

calculation.

Reference 4. The full-potential (FP) linear muffin-tin orbital (LMTO) calculation.

ct Reference

12. Using ultrasoft (US) pseudopotential and Vienna Ab-initio Simulation

Package (V ASP).

Reference 14. Using the linearized augmented plane wave and mixed-basis

pseudopotential methods.
f Reference 4.

Chapter3

3.4 Dislocation core structure revisited

3.4.1 Dislocation core energy variations with its polarization
In Section 3.2.4, we pointed out that the main difference between the asymmetric
and symmetric cores is their polarizations. The symmetric core does not have
polarization, while an asymmetric core has. In the previous ab initio calculations 3•12 , the
observed equilibrium screw dislocation has a symmetric core. However, our qEAM FF
predicted an equilibrium asymmetric screw dislocation core. To understand the difference
between the ab initio and the qEAM FF calculations in predicting dislocation core
structures, we calculated the dislocation core energy variation with its polarization both
using ab initio method and the qEAM FF.

We constructed a dislocation core with a particular polarization by translating the
central six columns of atoms along the [111] direction. The starting structure is the
asymmetric core (p=0.094 b) obtained using the qEAM FF. The final positions of those
six central columns of atoms are determined by scaling their displacement with the
dislocation polarization. To make the quantum calculations feasible, we choose the
parameters of simulation cells as X=3a[ll-2], Y=Sa[l-10] and Z=a[lll], which lead to
90 atoms per simulation cell in total. We further relax the quadrupole dislocation cells
while fixing the Z ([111]) positions of the central six columns of atoms for every
dislocation.

We also compared the resultant dislocation quadrupoles with the previous ab

initio calculation 3 ' 15 using the same simulation size. In that ab initio study, the

Chapter3

equilibrium dislocation with a symmetric core was obtained by minimization. We define
the averaged atomic deviation from the ab initio structure (denoted as "AI-structure") for
each simulation cell (contains 90 atoms) as Eq. (6).

/lt =

_i=_I_ _ __

(6)

here, trF is the X, y or z component of the position for an atom in the simulation cell
predicted by the qEAM FF, while t/' is the x, y or z component of the position for the
corresponding atom in the "AI-structure." The N stands for the number of the atoms in
the simulation cell and is 90 in our case. Note that the "AI-structure" has been scaled to
have the same equilibrium lattice parameter as those from the qEAM FF calculations.

G-E) X component

[3-£] Y component

Z compoueut
N= 90atom.'i

-----'---'------'-----'---'------'-----'---'------'-_J_-..L__

0.02

o.04

0.06

o.os

dislocation polarization (b)

0.1

o.'12

Chapter3

Figure 3-9. The averaged deviation of the atomic positions (x, y and z component)
between the 90-atom equilibrium dislocation quadrupoles from the qEAM FF and the
"AI-structure" from the ab initio (Ref. 3) calculations. In each equilibrium dislocation
quadrupole from the qEAM FF calculation, the dislocations have the same and specified
polarization. The "AI-structure" (from Ref. 3) contains four dislocations whose cores are
symmetric and zero polarized.

Figure 3-9 shows the averaged x, y and z component differences between the
qEAM FF predicted dislocation quadrupole with different polarizations and the "AIstructure." When polarization p = 0, the averaged atomic deviation between our resultant
structure and the "AI-structure" is small (~x = 0.004 A, ~y = 0.002 A, and ~z = 0.006 A).
It means that our resultant symmetric core structure is very close to the "AI-structure"

directly predicted using ab initio method. When the dislocation polarization increases,
the averaged atomic deviations of x and y components of position increase very little. ~x
changes from 0.004 A for the dislocations with symmetric core to 0.005 A for the
dislocations with fully polarized core; ~y changes from 0.002 A for the dislocations with
symmetric core to 0.003 A for the dislocations with fully polarized core. On the contrary,
the z component of atomic deviation increases from 0.006 A to 0.016 A (increases by
167%) when the polarization of the dislocations in the simulation cell changes from O b
(symmetric core) to 0.118 b (fully polarized core). These results are consistent with our
previous claim that the difference between the symmetric core and the asymmetric core
for a screw dislocation is mainly the dislocation polarization in the [111] (z) direction.

Chapter3

For each simulation cell, which contains a quadrupole of dislocations with the
particular polarization, we performed the one-energy evaluation using both the qEAM FF
and the DFf-LDA calculations 16 . Choosing the core energy of the symmetric core as the
reference point, we plot in Figure 3-10 the dislocation core energy variations with its
polarization.

A.b l11itio

Ta
;.,..
cJj

,.;

O.l

-73i 0.05

;.e;

qEAM FF

dislocation polarization (b)

Figure 3-10. The calculated dislocation core energy variations with its polarization using
the qEAM FF and ab initio methods.

Figure 3-10 shows that the dislocation core energy has a minimum (0.045 e Vlb
less than the symmetric core) at p = 0.094 b, leading to an equilibrium asymmetric
dislocation core structure in the qEAM FF calculations. While the ab initio computation
predicts that the symmetric or very close to symmetric dislocation core is energy

Chapter3

favorable, because the core energy is much higher when dislocation polarization is large.
Figure 3-10 gives us an explanation why the qEAM FF produces a different dislocation
core structure from the ab initio calculation. Though the ab initio methods are the most
rigorous way to evaluate the system energy, we still should bear in mind that the
employed 90-atom simulation cell is pretty small such that dislocation cores are
contacting with each other. The plot in Figure 3-10 might change if a much larger
simulation cell were used. At this moment, we consider our qEAM FF has its limitation
and need re-parameterizing to reproduce the ab initio results in Figure 3-10.

3.4.2 Force field re-parameterization
There are 19 tunable parameters in our EAM model force field. Our original force
field (denoted as "qEAMl ") was originally trained to optimally fit the following quantum
results 17 •
(1) The zero temperature equation of state (EOS) of Ta for bee, fee, and hep
crystal structures for pressures up to~ 500 GPa,
(2) the elastic constants,
(3) the volume relaxed vacancy formation energy also as a function of pressure,
(4) the equation of state for the A15 structure of Ta,
(5) the (100) surface energy in the bee Ta,
(6) the energies for shear twinning of the bee Ta.
To investigate the possibility whether the EAM model force field can predict the same
dislocation core structure as the ab initio calculations or not, we re-parameterized the

Chapter3

force field to fit the above quantum results as well as the quantum results of dislocation
core energy variations with its polarization.
During the force field re-parameterization process, we obtained another three
versions of the qEAM force fields and denoted them as "qEAM2," "qEAM3," and
"qEAM4." Table 3-3 gives the calculated lattice parameters a (A) and elastic constants
(Cn, C 12 and C44 in unit of GPa) for bee Ta model crystal.

Table 3-3. Experimental and theoretical values of lattice parameter a (A), elastic
constants Cn (0Pa), C12 (0Pa) and C44 (0Pa).

a (A)

Cn (0Pa)

C12 (0Pa)

C44 (0Pa)

qEAMl FF

3.32

273

138

69.6

qEAM2FF

3.35

255

148

60.2

qEAM3 FF

3.32

257

148

77.3

qEAM4FF

3.33

254

155

67.4

Ab Initio a

3.36

244

160

66.3

3.30

266

158

87.4

Exp.

a Reference 14. Using the linearized augmented plane wave and mixed-basis

pseudopotential methods. They are the inputs for force field fitting.
b Reference 13.

The results in Table 3-3 show that the force fields (qEAMl, qEAM2, qEAM3,
and qEAM4) all describe the bee Ta crystal well, however, Figure 3-11 shows that these
force fields lead to different behaviors of the dislocation core energy variation with its
polarization. The qEAMl FF predicted an equilibrium asymmetric dislocation core

Chapter3

structure and an energy minimum when the dislocation polarization p is 0.094 b. The
qEAM2, qEAM3 and qEAM4 force fields predict that the dislocation core energy would
increase with the increase of the dislocation polarization. The qEAM3 FF describes a
large rate of increase and best fit to the ab initio results, the qEAM4 FF describes a very
small increase rate, and the qEAM2 FF leads to an intermediate increase rate. Regardless
of the difference in the increase rate of the dislocation core energy with its polarization,
these three force fields all predict an equilibrium symmetric dislocation core for the
1/2a screw dislocation in Ta. The polarizations of the resultant dislocation core are
less than 1 x 10-4 b, which is at the same magnitude but smaller than the 7 x 10-4 b from
the MGPT FF calculation4 .

0.2

0 QM

0.15

,-_

..,

'l.)

qEM13

0.1

tJEAM2

:.a

Si: 0.()5

;.:..:

qEA:\14

qEA;\ll

-0.05

0.02

0.04

0.06

0.08

0.1

0.12

0.14

Polmization (b)

Figure 3-11. The calculated dislocation core energy variations with its polarization using
the qEAM force fields (qEAMl, qEAM2, qEAM3, and qEAM4) and ab initio methods.

Chapter3

We conclude here that the EAM model force field can predict the same symmetric
dislocation core structure as the ab initio method even though it does not completely
reproduce the quantum results (such as the qEAM4 FF). The above findings pose two
fundamental questions. (1) What is the reason for different force fields predict different
equilibrium dislocation core structures? (2) What is the affect of this difference to the
study of plasticity of bee Ta? In Section 3.4.3, we will report our preliminary results on
the first issue. The Question (2) is still under investigation. It should be noticed that the
reported calculation results in the following chapters are obtained in simulations using the
original qEAM FF unless specified.

3.4.3 Generalized stacking-fault energy (y) surface

The generalized stacking-fault energy (y) surface is an energy profile of two semiinfinite half crystals first displaced relative to each other by a vector v on a
crystallographic plane, then relaxed only in the direction perpendicular to the plane6 . The
y surfaces are the major input parameters to the well-developed generalized Peierls-

Nabarro model 18 , which is a continuum model for describing dislocation properties. The
generalized Peierls-Nabarro model informed with y surfaces calculated using the ab initio
electron theory has been used to study the <001> dislocations in bee metal Mo and Nb 19 ,
super-dislocations in Ni 3 Al 20 and NiA1 21 , and dislocations in fee metal Al 22 •23 . The y
surfaces are also considered very important for accurately modeling bee screw
dislocation behavior. Duesbery and Vitek 11 used the <111> cross section of the {110}
plane y surface in tantalum and molybdenum to explain the observed screw dislocation

Chapter3

core structure difference. They found for the l/2a screw dislocation a symmetric
core in Ta but an asymmetric core in Mo.

We calculated the and cross sections of the {ll2} y surfaces and
the cross section of the {ll0} y surface using the force fields qEAMl, qEAM2,
qEAM3, and qEAM4. These results are displayed in Figure 3-12 comparing with the ab

initio calculations4 . In our force field calculations, two parallel equivalent generalized
stacking-fault surfaces were introduced in a periodic bee Ta crystal cell. The distances
between two stacking-fault surfaces are 96 atomic planes for the {112} surface and 32
atomic planes for the { 110} surface, which are two times larger than the previous MGPT

FF calculations • After constructing the initial surfaces, we relaxed atoms in the direction
normal to the specified surface to the full convergence that the force on each atom is no
more than 3.5 x 10-4 eV/A. The ab initio result in Figure 3-12 is from Ref. 4. These
results are obtained by relaxing 12-plane supercell using the pseudopotential (PP)
techniques, then evaluating energies for the defined geometries using the full-potential
linear muffin-tin orbital (FP-LMTO) method.

Chapter3

0.25

qEA:\U:

L'.\ITO:
■ OD

0.2

.-.

~-<

>-'.::'.(l.15

{112}/<110>

;;,,
:,t;

....

: :l

0.1

{ll2}/
0.05

Normalized displacement (X

(a)

qEA:\12:

0.25

L:\ITO:

0 D

{112}/

°~.,, - -

0.05

, ,. _

.,,

.[J -

........ q_ 0

{110}/

.....

.._ ,._

- = - - - - ' - - - ' - - - ~ - ~ ~ - - - - ' - - - ' - - - ~ - ~ ~.• ,

0.2

0.4

0.6

Normalized displacement ex

(b)

0.8

Chapter]

qEA:\13:

0.25

L'.\ITO:
■ 0 □

Normalized displa,cmenl a

(c)

qEA:\14:

0.25

UHTO:
■ 0 □

';;!,r

0.2

-.::

,2,
>-..

;:f{J.15
·\)

-~:::

i"""1

0.1

{ll2}/

0.05

.-0:::

0 .

_,,/

,., -

~ - - ....

{H0}/<111>

.... 9.0

'..._

-"---,.__---'------'---'----'----'=,,,;.

0.2

0.4

0.6

Normalized displacement u.,

(d)

0.8

Chapter3

Figure 3-12. High-symmetry lines in the {112} and {110} y surfaces for bee Ta, as
calculated with the ab initio FP-LMTO electronic structure method and with the force
fields (a) qEAMl, (b) qEAM2, (c) qEAM3, and (d) qEAM4.

Figure 3-12(a) shows that the <111> cross section of both the {112} and the
{110} y surface from our original qEAMl FF quantitatively agree with the FP-LMTO
results but they surface {112)/<110> deviate severely from the ab initio results when the
normalized displacement is in the range of 0.3 ~a~ 0.7. On the other hand, Figure 312(b), (c) and (d) shows that they surfaces {112}/<111> and {110}/<111> calculated
using the qEAM2, qEAM3 and qEAM4 FF are lower than the ab initio results by almost
the same amount. More interestingly, the MGPT FF in Ref. 4 also made the same amount
of error. They surfaces {112 }1<110> calculated using the above three force fields agree
with the ab initio calculations except at a=0.5.

The qEAMl FF predicts an asymmetric core for screw dislocation, while the
qEAM2, qEAM3 and qEAM4 force fields predict symmetric cores for screw dislocation.
We find that the calculated y surfaces using these two groups of force fields show some
different features. This suggests that there is a possible relationship between the y
surfaces and the equilibrium dislocation core structure.

3.5 Conclusion

In this chapter, we first obtained the equilibrium dislocation core structure using
our qEAM FF. The qEAM FF leads to a polarized and asymmetric dislocation core,
spreading along three <112> directions in { 110} planes. We calculated the equilibrium

Chapter3

dislocation core energy both in atomistic model and continuum model. These two models
yield consistent results, 1.400 e Vlb in atomistic model and 1.404 eVlb in continuum
model. Furthermore, we obtained an insight to the dislocation core structure from an
energetic point of view. We are going to apply this insight to determine the dislocation
mobility in Chapter 4. By re-parameterizing the qEAM force field, we further show that
the EAM model force field is able to lead to an equilibrium symmetric dislocation core as
the ab initio calculations. Our calculated high-symmetry lines (<110> and <111>) in the
{112} and { 110} y surfaces for bee Ta using different groups (predicting asymmetric or
symmetric dislocation core) of force fields show different features. This suggests the
possible connection between the calculated y surfaces and the resultant dislocation core
structure.

Chapter]

3.6 References
1. A. M. Cuitino, L. Stainier, G. Wang, A. Strachan, T. <:;agm, W. A. Goddard, and M.
Ortiz, A Multiscale Approach for Modeling Crystalline Solids, J. Comput.-Aided
Mater. Design, to be published.W. Xu and J. A. Moriarty, Phys. Rev. B, 54, 6941
(1996).
2. W. Xu and J. A. Moriarty, Phys. Rev. B, 54, 6941 (1996).
3. S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett., 84, 1499 (2000).
4. L. H. Yang, P. Soderlind, and J. A. Moriarty, Philos. Mag. A, 81, 1355 (2001).
5. J. P. Hirth and J. Lothe, Theory of dislocations (Wiley, New York, 1982), p.465.
6. V. Vitek, Cryst. Lattice Defects, 5, 1 (1974).
7. A. Seeger and C. Wuthrich, Nuovo Cimento, 33B, 38 (1976).
8. W. Xu and J. A. Moriarty, Comput. Mat. Sci., 9, 348 (1998).
9. J. A. Moriarty, W. Xu, P. Soderlind, J. Belak, L. H. Yang and J. Zhu, J. Eng. Mater.
Tech., 121, 120 (1999).
10. W. Sigle, Philos. Mag. A, 79, 1009 (1999).
11. M. S. Duesbery and V. Vitek, Acta Mater., 46, 1481 (1998).
12. C. Woodward and S. I. Rao, Philos. Mag. A, 81, 1305 (2001).
13. Single Crystal Elastic Constants and Calculated Aggregate Properties: A Handbook,
edited by G. Simmons and H. Wang (MIT Press, Cambridge, MA, 1971).
14. 0. Gi.ilseren and R. E. Cohen, Phys. Rev. B, 65, 064103 (2002).
15. S. Ismail-Beigi, D. E. Segall and T. A. Arias, private communication.

Chapter 3

16. Using the program seqQuest and LDA pseudopotential for Ta developed by P. A.
Schultz.
17. A. Strachan, T. <;agm, 0. Gtilseren, S. Mukherjee, R. E. Cohen, and W. A. Goddard,
First Principles Force Field for Metallic Tantalum, Phys. Rev. B, submitted.

18. G. Schoeck, Philos. Mag. A, 69, 1085 (1994).
19. G. Schoeck, Philos. Mag. Lett., 76, 15 (1997).
20. G. Schoeck, S. Kohlhammer and M. Fahnle, Philos. Mag. Lett., 79, 849 (1999).
21. G. Schoeck, Acta Mater., 49, 1179 (2001).
22. G. Lu, N. Kioussis, V. V. Bulatov and E. Kaxiras, Phys. Rev. B, 62, 3099 (2000).
23. G. Schoeck, Philos. Mag. A, 81, 1161 (2001).

Chapter4

Chapter 4 Peierls energy barrier and Peierls stress of 1/2a screw
dislocation
4.1 Overview
To understand plasticity of crystals, it is critical to understand the mobility of
dislocations and the role of the dislocation core in the slip process. The Peierls-Nabarro
model provides an analytical strategy to compute the required stress to move a
dislocation in an otherwise perfect crystal from the misfit energy. This model suggests
that during the translation of a dislocation, the entire variation in the potential is
associated with the changes in the shear misfit energy between two half-crystals and with
no variation in the elastic field of the dislocation 1 . The Peierls-Nabarro model applies
even though the dislocation core reconfigures during the motion 2 • In previous atomistic
simulations 3' 4 , the Peierls stress was determined by increasing the applied shear stress
incrementally and fully relaxing the simulation cell (containing a dislocation) until the
dislocation glides. In the following, we present an alternative approach to obtain Peierls
energy barrier and Peierls stress directly from the analysis of a moving dislocation, which
we achieve by simulating a dislocation dipole migration process at extremely low
temperature (0.001 K).

As pointed out in Ref. 5, the twinning and anti-twinning slip asymmetry (details
see Section 4.4) of shear on { 112} planes is an intrinsic factor in the observed violation
of the Schmid law for plastic behavior of bee metals. Owing to this asymmetry, the
required shear stresses along [111] and [-1-1-1] in {112} planes for the same dislocation
to glide are not equivalent. By measuring the dislocation core energy as a continuous

Chapter4

function of position in the lattice, we observe this twinning and anti-twinning slip
asymmetry of the Peierls energy surface and dislocation motion trajectory. In this
chapter, we also report our studies on dislocation motion at finite temperatures and via
kink pair mechanism.

4.2 Dislocation dipole migration and annihilation process

We construct a dislocation dipole (two dislocations with opposite Burgers
vectors) using elastic theory from the perfect crystal with periodic boundary conditions.
The simulation cell contains 5670 atoms with lattice vectors X=9a[ll-2], Y=15a[-110]
and Z=7a/2[111]. In the simulation cell, the dislocations with b=a/2[111] and b=a/2[-1-11] are positioned at (l/2X, 3/4Y) and (l/2X, 1/4Y) in the (111) plane, respectively. This
is denoted as the [1-10] dislocation dipole. Keeping the lattice parameters fixed to the
perfect crystal values, the introduction of this [1-10] dislocation dipole causes stresses of
O'xz = -1080 MPa, O'xx = 410 MPa, O'yy = 530 MPa, O'zz = 250 MPa and O'xy = 0, O'yz = 0. The

large xz shear stress is due to the misfit of atomic positions in the (1-12) cell boundary
(see Appendix). Since the non-glide stresses could have effect on the computed Peierls
stress 3,4· 6 , we first relaxed the stress of the simulation cell using NPT MD simulations
(with the Rahman-Parrinello barostat7 and the Hoover thermostat 8 at T = 0.001 K) to
reach a zero stress [1-10] dislocation dipole. The final lattice parameters are given in
Table 4-1. The volume increase for a dislocation per Burgers vector from the bulk lattice
cell to the zero-stress cell is 4.29 A3 .

Chapter4

Table 4-1. The lattice parameters and volumes for simulation cells (containing the [1-10]
dislocation dipole) with zero stress, bulk lattice parameters and different pure shear
stresses. a denotes the angle between axes Y and Z; ~ denotes the angle between axes Z
and X; while y denotes the angle between axes X and Y.

IXI

IYI

IZI

Volume
(A3)

(A)

CYxz = 500 MPa

73.38

70.54

20.13

90.00

91.59

90.00

104178

CYxz = 300 MPa

73.37

70.54

20.13

90.00

91.39

90.00

104172

axz=0MPa

73.36

70.55

20.13

90.00

91.12

90.00

104170

CYxz = -300 MPa

73.35

70.56

20.13

90.00

90.84

90.00

104168

axz = -500 MPa

73.34

70.56

20.13

90.00

90.67

90.00

104167

CYxz= -1100 MPa

73.34

70.57

20.13

90.00

90.13

90.00

104174

bulk lattice parameters

73.27

70.51

20.15

90.00

104110

In order to study dislocation migration, we applied an external shear stress to our
simulation cell (in the twinning and anti-twinning directions). In a [1-10] dislocation
dipole, both dislocations are sheared in the (1-12) plane in the twinning sense under the
shear stress along the [111] direction (CTxz > 0) and in the anti-twinning sense when the
shear stress is in the [-1-1-1] direction (O"xz < 0). In these two cases (O"xz > 0 and O"xz < 0),
we start from zero stress and then increase (twinning) or decrease (anti-twinning) the
applied shear stress O"xz on the simulation cell in steps of 100 MPa until the dislocations
begin to move. For each stress state, we performed 10 ps of NPT MD simulation
followed by 25 ps of NVT MD simulation at 0.001 K. We find that the dislocation dipole

Chapter4

starts to move under twinning shear of O'xz = 500 MPa and anti-twinning shear of O'xz = 1100 MPa.

Once the dislocation dipole starts to move under the applied shear stress, we
continued the NVT MD simulation up to 125 ps. In the course of simulation, dislocations
move continuously until the annihilation occurs. The solid lines of Figure 4-l(a)
(twinning) and Figure 4-l(b) (anti-twinning) show the time evolution of the total strain
energy (the sum of the atomic strain energies calculated using Eq. (2) in Chapter 3 with
reference to perfect crystal) per dislocation per Burgers vector during the dislocation
migration and annihilation process. The total strain energy decreases as the dislocations
approach each other. The rapid drop of the total strain energy at the end indicates
dislocation annihilation. Figure 4-l(a) for the twinning shear (O'xz = 500 MPa) shows a
residual total strain energy of 0.2 eV (per dislocation per b) after dislocation annihilation.
This is because the initially set cell parameters are different from the final dislocation free
crystal cell lattice parameters (cf Table 4-1 ). However, there is little residual strain
energy in Figure 4-l(b) for the anti-twinning shear (O'xz= -ll00 MPa), because the lattice
parameters of the simulation cell under O'xz = -1100 MPa are very close to those of perfect
crystal as shown in Table 4-1. This is reasonable because building the [1-10] dislocation
dipole into crystal using the lattice parameter of the perfect crystal leads to a shear stress
of O'xz = -1080 MPa.

Chapter4

11 -1 0) dislocation di1>olc under twinning shear
- - - - total strain cnerg_-v
. • • . • . • ela~tic energy

o:xz = 500 MPa

__ -

- core cnerg_v

...
o------'---'--_,____ ___.__ __,_ _

J.....__

_j__----='-'-""---.......~ -

1.25

100

time (ps)

(a)

[l -1 OJ dislocation dipole under antitwlnnlng shear
- - - - total strnin encrg_y

crxz = -1100 MPa

• • • • • • • elm;tie energy

- core energy

75
time (ps)

(b)

100

125

Chapter4

Figure 4-1. The variation of the total strain energy, elastic energy, and core energy with
time in the NVT MD simulations at T=0.001 K. These simulations simulate the migration
and annihilation of the [-11 0] dislocation dipole under the smallest shear stress required
for dislocation migration. (a) Twinning shear ( Dxz = 500 MPa) and (b) Anti-twinning
shear (O"xz = -1100 MPa). There are 5670 atoms in the periodic simulation cell. The states
A and C correspond to the minimum core energy configurations, while state B
corresponds to the maximum core energy configuration as shown in Figure 4-3. The
detailed structures for these states are shown in Figure 4-2 using the corresponding DD
maps.

The dynamics in Figure 4-1 show that the total strain energies exhibit bumps on
top of the generally monotonic decrease as the dislocations in the dipole migration. To
understand the origin of these bumps, Figure 4-2 shows the DD maps of the dislocation
dipole for the points labeled as A, B, and C in Figure 4-1. Panels (A) and (C) in Figure 42(a) and Figure 4-2(b) show that the valleys of the energy bumps have configurations in
which the dislocations are in equilibrium positions. In contrast, panel (B) shows that the
peak of the energy bump corresponds to a configuration in which the dislocation is
halfway between two equilibrium positions. Thus, the bumps in total strain energy relate
to the dislocation motion through a periodic Peierls energy barrier of lattice resistance.
Figure 4-2 also shows that during each step the dislocation moves by a/3<112> on { 110}
planes regardless of the sense (twinning or anti-twinning) of shear. This leads to a zigzag
path for dislocation motion as shown by the dotted lines in panel (D).

Chapter4

500 MPa)

(a) Twinning (O'xz

.. ◄ o•• ◄,. a •
' • • o
o • '• "•o

'"'. Q "" ... ~ ,·,
\ . : •• 0

• I

O►

• '<

¼0- ►

,.<

lic ► O•~O

(D)

[-110]

[ 111]

(b) Anti twinning (O'xz

------►

[11 -2]

-1100 MPa)
(D)

"o•'?o•"o
•~o• ■

-c,i,o•

.011~''·.~.:.
"o ,•

"tf • * o

►o

' "

.... {i,--Q ►•
... ... ,,

Q-

Iii ◄
• •

,._

• ► <1> ► 0>•>"'-

\ ' ,. ► o,., •
• ~ ,. ... a ....

...+

4 • o,,.f~>-:0"' • C..,

[-11 0]

•·..•..

.....•····

t 0

..··.
..•'

-:::...

..·
..·..
·::....
-..o.,,,o

....

[ 111]

------►

[11 -2]

Chapter4

Figure 4-2. DD maps of the states of system during the dynamical process of dislocation
dipole migration and annihilation under (a) twinning and (b) anti-twinning shear. Panels
(A) and (C) show the dislocation dipoles at equilibrated states, while (B) show the
dislocation dipoles translating halfway between two equilibrium states. These maps show
only the central region of the simulation cell containing 5670 atoms. Panels (D) show the
[111] projections of dislocation dipole. The dotted line in (D) plots the dislocation slip,
which is in a zigzag style along <112> directions in {110} planes, in the simulation. The
arrow beside dislocations (drawn as plus sign and minus sign) in panel (D) indicates the
direction of the Peach-Koehler force for that dislocation introduced by the applied shear
stress. The directions of the shear stresses O"xz are represented by the plus sign ([111]) and
minus sign ([-1-1-1]) in the circles in panel (D).

In Section 3.3.1 we defined the dislocation core as the twelve atoms with highest
strain energy per Burgers vector for the equilibrium dislocation. We now apply this
definition to any configuration of a dislocation during its motion. In this way, the total
strain energy of our system can be partitioned into two parts: core energy and elastic
energy. The dashed lines in Figure 4-1 show the time evolution of the core energy and
the dotted lines show the similar curve for the elastic energy. The core energy is rather
constant throughout the dislocation migration process, showing bumps at the same places
as the total strain energy. The core energy rapidly drops to zero as the dislocation pair is
annihilated. This clear-cut definition of elastic energy leads to desirable smooth
monotonic decrease with no bumps as the two dislocations move towards each other.
These results support the validity of our definition of dislocation core energy.

Chapter4

4.3 Peierls energy barrier and Peierls stress from dislocation dipole migration

1.54

[ I -1 OJ dislocation dipole under twinning shear

1.52

~ l.48

a/2[-l -1 - ll dislocation
a/2[1 l ll di'ilocation
cosine function fitting

""
f::""

§ 1.44
::::

i:J

~ 1.42

1.38

136•0--~-2~-~--4~-~-~6--'---~s-~----w
translation distance(,\)

(a)

1.54

(1 -1 OJ dislocation dipole under antitwinnlng shear

o a/2[-l -1 -lJ dislocation
l..52

a/2[1 l IJ dislocation
cosine function fitting

1.4

1.38

l.36---~-~--'---'------'L._--'----'---'--_L---1...___
10

translation distance (i\)

(b)

Chapter4

Figure 4-3. Dislocation core energy as a function of distance traveled by the dislocation
under (a) twinning and (b) anti-twinning shear. The solid line shows the cosine function
[Eq. (1)] fit to the atomistic data. Table 4-2 gives the parameters from optimal fitting. The
states, denoted as A, B and C, correspond to those shown in Figure 4-2.

We define the dislocation position as the strain-energy weighted geometric center
of the 12 atoms forming the dislocation core per Burgers vector. Figure 4-3(a) shows the
variation of dislocation core energy with the dislocation translation distance for the case
of twinning shear and Figure 4-3(b) shows the similar plots for the case of anti-twinning
shear. Both core energy variation curves fit well the following cosine function in Eq. (1).
2n:x
E (x) = _ P [1-cos(-+ m)] + k · x+ Ee,
'f'

(1)

Here xis the distance traveled by the dislocation and Ec(x) is the dislocation core
energy at translation position x. The parameter Ep is the Peierls energy barrier, L is the
translation distance for a single dislocation jump, and Ee is the dislocation core energy at
its equilibrium position. A phase shift (fJ and linear term k·x are also introduced to better
describe our data. Table 4-2 gives the fitting parameters for dislocation motion under
twinning and anti-twinning shears. The Peierls energy barriers are determined to be
Ep(twinning)=0.032 eV/b and Ep(anti-twinning)=0.068 eV/b. The anti-twinning to
twinning ratio of Peierls energy barrier is Ep(anti-twinning)/Ep(twinning) = 2.125.
Table 4-2. The parameters obtained from fitting the dislocation (lb long) core energy to a
cosine function [Eq. (1)] of its translation distance under twinning and anti-twinning
shears. The parameters include the Peierls energy barrier Ep in eV, periodic translation

Chapter4

distance Lin A, dislocation core energy E, in eV, phase shift term.

Ep (eV)

L (A) a

E, (eV) b

k (eV/A)

Twinning

0.032

2.48

1.414

0.83

0.003

Anti twinning

0.068

2.90

1.401

1.93

0.003

Compare with la/3<112>1 = 2.72 A in perfect crystal.

b Compare with the equilibrium dislocation core energy E, (eq.) = 1.400 eV.

In the Peierls-Nabarro model, the stress z(x) felt by the dislocation during its

motion is the derivative of the core energy with respect to the distance it traveled 1.

_ 1 dEc (x)
b 2 dx

( )TX

(2)

Substituting only the cosine term in Eq. (1) to Eq. (2) leads to the Peierls stress
[the maximum stress from Eq. (2)] in Eq. (3).

1 n ·E
b2 L

(3)

Using the Ep and L obtained above (Table 4-2) and the Burgers vector b=2.88 A,
we determine Peierls stresses of
;(twinning)=790 MPa or ;(twinning)/µ= 0.013,
;(anti-twinning)=1430 MPa or ;(anti-twinning)/µ= 0.024,
;(anti-twinning)/;(twinning) = 1.80.
Here, the calculated shear modulus for the perfect crystal of µ=62.3 GPa 9 .

Chapter4

An alternative approach for calculating the Peierls stress [also called the critical
resolved shear stress (CRSS)] is to shear an infinite cylinder containing a dislocation.
Using the same qEAM FF the calculated CRSS=740 MPa ' 0 is in good agreement with
the 790 MPa derived above for twinning shear with periodic boundary conditions.
Table 4-3. The computed Peierls stresses ; in unit of MPa for twinning and antitwinning shears for l/2a screw dislocation in Ta.

F-S a

MGPTb

qEAM

(Ito et al.)

(Yang et al.)

(present work)

0.0007

0.09

; (twinning)

4120

600

790

T,z (anti-twinning)

14800

1380

1430

Anti-twinning/ Twinning ratio

3.59

2.29

1.80

Force Fields

Dislocation Polarization (b)

Reference 4. The reported ; is 0.05 C44 for twinning shear and 0.18 C44 for anti-

twinning shear. To calculate ; in MPa, we used C44 = 82.4 GPa from Ref. 5.

b Reference 3. The reported ; is 0.0096 G for twinning shear and 0.022 G for antitwinning shear. The shear modulus G is 62.5 GPa.

Table 4-3 compares our results for Peierls stresses with previous calculations
using other force fields. The recent calculations using the MGPT potential 3 lead to 600
MPa for twinning (24% less than ours) and 1380 MPa for anti-twinning (3% less).
Calculations using the simple Finnis-Sinclair (F-S) potentia1 4 lead to a nonpolarized
dislocation and much different values (4120 MPa for twinning and 14800 MPa for anti-

Chapter4

twinning). Thus, even though qEAM FF leads to a larger dislocation polarization than
MGPT FF, the two force fields lead to a similar description of dislocation mobility.

Our best estimate of the Peierls stress (790 MPa) at 0.001 K is still larger by a
factor of two than the best extrapolation to OK (-300 MPa 11 ) from experiments at finite
temperatures. A possible reason for this disagreement is that dynamic kink-like processes
might lower the Peierls stress at the experimental temperatures (> 73 K). Our simulation
at 0.001 K would not include such processes (due to the short length of the periodic cell
in the [111] direction). Section 4.5 considers simulations of the dislocation motion at
finite temperatures, where we find Peierls stresses in reasonable agreement with
experiment.

4.4 Twinning/Anti-twinning asymmetry

In bee crystal, there is the inherent twinning and anti-twinning asymmetry in
{ 112} planes 5 ' 12 • The stacking sequence of { 112} planes has a six-layer repetition
... ABCDEF .... A displacement a/6[111] on (-1-12) plane produces a stacking sequence of
... ABCDCDEFAB ... , which corresponds to a monolayer twin and it is natural to regard
this as a possible single layer fault (twinning fault). However, the displacement in the
opposite sense (-a/6[111]) would produce a stacking ... ABCDABCDEF ... , which is
different from the former one and does not correspond to a monolayer twin and thus the
fault created is different from the above "twinning" fault. It is called the "anti-twinning"
fault.

Chapter4

Our results show clearly the twinning and anti-twinning asymmetry of shear for
1/2a screw dislocation in bee Ta. We find an anti-twinning/twinning ratio of 2.125
for Peierls energy barriers and 1.80 for Peierls stresses. Table 4-2 also shows that the
derived dislocation core energy £,(twinning) = 1.414 eV/b is 1% higher than E,(antitwinning) = 1.401 eV/b. Both agree quite well with the dislocation core energy Ec(eq.) =
1.400 eV/b obtained by summing the atomic strain energies for the 12 atoms in the
equilibrium dislocation core. The difference between the periodic translation distance for
twinning shear (2.48 A) and anti-twinning shear (2.90 A) suggests that dislocations move
differently in these two cases.

Figure 4-4 shows the trajectories for dislocations with b=a/2[-1-1-1] [Figure 44(a)] or b=a/2[111] [Figure 4-4(b )] under twinning and anti-twinning shears. The
dislocation position in the figures is determined as the strain-energy-weighted geometric
center of the 12 atoms constituting dislocation core. The origin of the plot is the initial
position for the dislocation. Figure 4-4 shows that dislocations with b=a/2[-1-1-1] and
b=a/2[111] behave similarly under the same sense (twinning or anti-twinning) of shear,
while a dislocation moves along completely different trajectories under different senses
of shear. Under anti-twinning shear, the dislocation moves along a path at an angle of
29.5° with the [-110] direction. This angle is close to the 30° for the observed slip system
(<112> directions on { 110} planes) from DD maps. Because the dislocation trajectory is
not a straight line, the periodic translation distance 2.90 A of this path is larger than
la/3<112>1=2.72 A. However, for twinning shear the path of the dislocation makes an
angle of only 8.5° with the [-110] direction, leading to a shorter periodic translation
distance (2.48 A).

Chapter4

b = 1/2a[-1-1-1]

Antitwinning

·='E

-0

(1.=

(0-11)

}-(-110)

29.5°

()
~ = 8.5 0

-2

(-101)

Twinning

-6

-8

• JO

-8

-6

-4

-2

11 l -2] direction (..\)

(a)

b = 1/2a[l 11]

.~_,

-0

Twinning

8.5 0

= 29.5°

(-l O I)

}-(-IIOJ

()

-2

(0-1 H

An lit winning

-4
-6

-8
-101---'--'---'-'--'----'----'--'---'-___,_-'---'--'---'--'---'--'---'-"-a
-8
-6
-4
-2
IO
12

I1 l -2) direction (i\.)

(b)

Chapter4

Figure 4-4. The <111> projection of the motion trajectory for a dislocation with (a)
b=a/2[-1-1-1] and (b) b=a/2[111] under twinning and anti-twinning shears. The origin
represents the position of the initial equilibrium dislocation. The schematic map on the
right shows the crystal geometry and the twinning or anti-twinning direction of shears.
The path which dislocation follows under anti-twinning shear makes an angle of a=29.5°
with the [-110] direction while the path which dislocation follows under twinning shear
makes an angle of 13=8.5° with the [-110] direction.

Table 4-3 shows that the qEAM FF calculations (screw dislocation with a
polarization of 0.09 b) lead to ;(anti-twinning)/;(twinning) = 1.80, whereas the MGPT
FF calculations 3 (screw dislocation only slightly polarized 0.0007 b) lead to ;(antitwinning)/ ;(twinning) = 2.29 and the F-S FF calculations4 (screw dislocation with a fully
isotropic core, zero polarization) lead to ;(anti-twinning)/;(twinning) = 3.59. Thus, the
;(anti-twinning)/;(twinning) is larger than one and seems to increase as the polarization
decreases.

4.5 Dislocation motion at finite temperatures

The previous section used the Peierls-Nabarro model to analyze the dislocation
motion at T=0.001 K. However, at high temperatures, the thermal energy fluctuations are
too large for this atomistic-energy-based analysis to be useful. Hence, we now use
constant temperature and pressure (NPT) MD simulations to study the dynamic processes
of dislocation motion at finite temperatures. This study employed a larger simulation cell

Chapter4

containing 22,680 atoms and performed the simulations at temperatures (20K, 50K,
lO0K, and 300K) and under zero pressure.
For the cases of T = 20K and T = 50K, we find that dislocations always move in
<112> directions on (110) planes (the same slip system found in the simulations at
T=0.00lK). However, dislocations do not always move towards each other, sometimes
they take steps perpendicular to the dipole direction.

At 100K, the dislocations still move in the same slip system but the thermal
energy is large enough for the polarizations to occasionally change the sign without
jumping. Such processes were not observed in lower temperature simulations (T=0.00lK,
20K and 50K). Figure 4-5(a)-(d) shows the DD maps for the motion of one dislocation at
T=lO0 Kand time t=l.6 ps, 2.4 ps, 3.2 ps, and 4.4 ps. We observe that from t=l.6 ps to
t=2.4 ps the dislocation moves in the [-1-12] direction (left) [Figures 4-5(a) and Figure 45(b)] and changes its polarization after the hop/jump. From t=2.4 ps to t=3.2 ps [Figures
4-5(b) and Figure 4-5(c)], the dislocation changes its polarization but it stays at the same
position. Finally, at t=4.4 ps [Figure 4-5(d)], the dislocation moves forward to the next
equilibrium position along the [-12-1] direction and again changes polarization in the
jump. Figure 4-5(e) summarizes schematically the motion of this dislocation from t=l.6
ps to 4.4 ps.

Chapter4

•-o •

•-•

f \

o-• ◄ •-o

e-o

e-e ► O-e

o-•

•' o •

- ◄ ,!;

411

O••

411

••

411•0

••

e ► O

lll

• ►

0 ► 0

• •

o•

-•··

••0

• 0

" "
• ' o-•-•

O+e+lll•O

•••◄ o•• ◄•

, \ I

411

• • '

••// .'o'o/t1:~• ■

••

0 ••••O

~(_d_)_t_=_4_._4_~1: : : • ~ ...•0 • •

•••
• • ~,0 ,,;
•/,.
• •' o- "-o-•
' '
• • o •-•-o••-•-o • •

•• ,~,10,- •, •

•-o • • o • • o • • o
• • o • • o-• • o • •

o •

• •

0 - •

'•

(e) Dislocation motion at T=lOO K

(d) t

4.4 ps

1.6 ps

2.4 ps

Chapter4

Figure 4-5. Dislocation dipole dynamics process for a system containing 22,600 atoms at
T=lO0 K simulated with NPT MD. Snapshots at different times [(a) t=l.6 ps; (b) t=2.4
ps; (c) t=3.2 ps and (d) t=4.4 ps] are shown using DD maps. Only the 30x30 A2 region of
interest is shown. In the figure, the dislocation moves one step towards the right from 1.6
ps to 2.4 ps and changes polarization while staying in the same equilibrium position at
3.2 ps. This change of polarization allows the dislocation to move to the position shown
in figure (d). (e), The schematic representation of the dislocation motions from figure (a)
to (d).

The simulation at T=300 K also sometimes shows dislocation polarization
changes without a dislocation jump. Thus, Figures 4-6 (a)-(f) show the DD maps for a
dislocation at T=300 K. We observe that the dislocation changes its polarization from t=2
ps [Figure 4-6(d)] to t=2.5 ps [Figure 4-6(e)]. The process for this motion is sketched in
Figure 4-6(g).

Chapter4

(a)t=0.5ps
, '

• •

O ► e

• •

•-o

.._

• • ' 0 ' I· - .

• • , • •. . ' I " o • •
, •
• •
• •
• • a••
• •
• • •
• •
• •,
• • •-•
• • • • • • •
• 0 • •
• • • •

•••◄ o-•

•◄ O••

O ◄ e

....

--~ • • •

•◄ O

o-• ◄•

0 -

• • 0 • ' • ,o-•I •-o • •
• 0 •
, -o' •
• • 0 ,• , ,• ,o'
\1,,, •
• 0 • • 0 • •'•,• •, ' ,• 0
• • 0 '•' • o-•••, 0 ..,• " ,•
• 0
• 0
• • o-• ,• • •
•-o
• • 0 • • 0
• • o-• • 0 • • 0 • •
• o-• ◄• 0 • • 0 • •

•·•' ~
I \ ,o-•

.".. . o-•
t. ••o-•-•

(b)t=0.7ps

c:::-,

(c)t=l.2ps_,r.o.•o••o·• ►•

_______
L.

' . "o •
• •
,,,, '

111•0

• I

• ◄ 4l 1 -

o-•

/0\ : , 1 f/lJ

,. I

o-•
,. , - •

o-• ◄•◄ o••

•••◄ O

e-ill

o •••

o-• ◄•

O ◄•

o-•

• •

•-O

o •·•

o •

r. o,: •

(f) t = 2.6 ps

o-• ◄ I •-o

•-•-o • •

• ◄ O

-•

O;z_,•

••
► I
0•

o-• ►•

41ll~0

0 : • 0 •• 0 •

ill

' ,
•-o

•◄
O ◄ e

•-•

,•

O ◄•

•••

' • ...\ •

•◄-0

O ◄•◄•

•-•-0 •◄•

•-o ◄•

•·• ◄ C

Chapter4

(g)

Dislocation motion at T=300 K

(f) t = 2.6 ps

(e)t=2.5ps
(d)t=2ps

Figure 4-6. Dislocation dipole dynamics process for a system containing 22,600 atoms at
T=300 K simulated with NPT MD. Snapshots at different times [(a) t=0.5 ps; (b) t=0.7
ps; (c) t=l.2 ps; (d) t=2.0 ps; (e) t=2.5 ps and (f) t=2.6 ps] are shown using DD maps.
Only the 30x30 A2 region of interest is shown. A process of changing polarization of
dislocation can be seen in figure (d) and (e). (g), The schematic representation of the
dislocation motions from figure (a) to (f).
For the finite temperature simulations of dislocation dipole dynamics, we
compute the hopping rate (Y]) of dislocations as a function of temperature. We calculated
the hopping time as the average duration for the first 8 jumps for each dislocation and
then took the reciprocal to obtain the hopping rate. Figure 4-7 shows that the hopping rate
follows an Arrhenius behavior with temperature, yielding activation energy of 0.0053
eV/b for dislocation hopping. This activation energy is 6 times lower than the Peierls
energy barrier 0.032 eV/b obtained at 0.001 K for twinning motion.

Chapter4

lO ~~~-~~--.----r-------.--.-----.-----,--,--------,

qEAM Ta

Dislocation dipok
T=300 K

N=22,680

·-e l
....
T=IOO K
·o..( ;fl

c..

,.....
' l)
__,

~JJ

c..
,..
_,

0.01

0.02

0.03

0.04

0.05

0.06

I/temperature (l/K)
Figure 4-7. The Arrhenius plot of logarithm of dislocation hopping rate with reciprocal of
simulation temperature. The dislocation hopping time is defined to be the average
dislocation jumping duration for the first 8 jumps. The dislocation dipole containing
N=22,600 atoms was simulated with NPT MD at 20 K, 50 K, 100 K and 300 K,
respectively. This leads to activation energy for dislocation hopping of 0.0053 eV.

To estimate the Peierls stress at high temperatures, we assume that the Peierls
stress is proportional to the height of the energy barrier. This leads to an average Peierls
stress of Ip = 700x(0.0053/0.031) = 120 MPa for temperature in the range of 20K to
300K. This estimated average stress from the dynamical simulations agrees well with the
empirical average flow stress 11 (~110 MPa) in the same temperature region (170 MPa at
73 Kand 50 MPa at 273 K).

Chapter4

The dramatic drop of the activation energy for dislocation motion from 0.032
eV/b at 0.00IK to 0.0053 eV/b for the range of 20K to 300K implies a change in the
nature of the dislocation dynamics between these temperature regions. At 0.00IK, we
find that the dislocation moves collectively as a whole with every part of the dislocation
overcoming the same Peierls energy barrier. However, at finite temperatures, thermal
fluctuations cause different segments of the dislocation to move in a less correlated way
(but without creating a dislocation kink).

To study the effect of temperature on the dislocation topology, we heated the
dislocation quadrupole (5670 atoms) at temperature increments of 25K. At each
simulation temperature, we carried out 10 ps NPT MD (the pressure is 0 GPa) to allow
the volume to change and followed by 25 ps of NVT MD. At this heating rate, we find
that the dislocations fluctuate thermally about their centers but do not migrate.

analyzed the dislocation polarization at different temperatures and found that the average
polarization of the dislocation does not change with temperature. However, as shown in
Figure 4-8(a), the standard deviation of the dislocation polarization fluctuation increases
with temperature and is about 8 times larger at 300K than at IK. Figure 4-8(b) shows the
correlation coefficients for the polarization fluctuation between two nearest neighboring
one Burgers vector long dislocation pieces at different temperatures. These results
indicate that the correlation coefficient between the neighboring pieces of a dislocation
decreases with increasing temperature. Apparently, this decreased correlation may
decrease the effective energy barrier for dislocation motion.

Chapter4

dislocation quadrupole
N=5670

0.9
0.8

0.7
-.::;

-...

e 0.6
::::

.'2
C:

":... 0.5

c,;

1.)

,_,

0.3
0.2
0.1

100

150

200

250

Temperature (K)

(a)

minimi1..ation at OK

dislocation <1uadrupolc
N=5670

0.9

0.7

100

150

Ternperatme (K)

(b)

200

250

300

Chapter4

Figure 4-8. (a) The average polarization of the dislocation core as a function of
temperature. The error bars in the figure indicate the standard deviation of the dislocation
polarization. (b) The correlation coefficients of the dislocation core polarization between
the first nearest neighboring dislocation segments as a function of temperature. The data
at OK is derived from a minimization simulation and the others are computed by
analyzing the last 15 ps simulation trajectory from a total 25 ps TVN MD simulation.

4.6 Dislocation motion by nucleating kinks

In previous sections, we studied the dislocation motion in the cases that it moves
as a whole line. We are also interested not only in the case that the dislocation moves via
kink pair mechanism, i.e., nucleation of kink pairs and propagation of kinks along the
dislocation. The process of kink pair nucleation is hard to study in atomistic scale
because the relatively short lengths of dislocations (~ hundreds of Burgers vectors) and
short times (~hundreds of picoseconds) make the kink pair nucleation along a dislocation
very unlikely to happen during the course of a MD simulation. To remedy it, we studied
the process of dislocation migration and annihilation for a [1-10] oriented dipole of screw
dislocations via MD at T=O.OOlK and provided a nucleation center for the kink pair by
introducing a vacancy in the path of a dislocation (denote as Dv)- The other dislocation
(Df) sees a defect-free environment. We used a relatively long simulation cell (N=56,700
atoms) whose lattice parameters are X=9a[l-12] (73.356 A), Y=15a[l-10] (70.548 A)
and Z=70a/2[111] (201.5 A). We find that dislocations and vacancies attract each other.
For Ta, the vacancy formation energy in the core of a screw dislocation is 2.45 eV, while
the vacancy formation energy in the perfect bee crystal is 2.95 eV.

Chapter4

l.52r--a----,-r-,--,---,---,--,-,--,-,---,--,---,----,---,-----,------,----,--,----,----,---.----,

Ta
dislocation with vacancy (D )

kink formation

;> lA-8
3,
;:,-,

ff lA-6

0,)

lA-4
lA-2
lA-.__.__,___,_~~___.__,__.__~~~~__,_~~___.____.__,__.__~~'--'

time (ps)

Figure 4-9. Dislocation core energy as a function of time for dislocation dipole migration.
The full line shows the dislocation with a vacancy in its path will nucleate a kink pair to
reduce the activation energy for a jump.

From our MD simulations of dipole migration, we find that the vacancy helps the
dislocation nucleate a kink pair that propagates making the dislocation advance. In Figure
4-9, we show the core energy per Burgers vector as a function of time for two
dislocations (Dv and Df). The core energy is defined as the strain energy of the 12 atoms
with higher atomic strain energy per dislocation, per Burgers vector. The dislocation
without the vacancy (Df) shows the same motion behavior we found in smaller (7b long)

cells 1 , shown as the dashed line in Figure 4-9. It moves as a rigid straight line with an
activation barrier of 0.07 eV/b. On the other hand, the dislocation with the vacancy in its

Chapter4

way (Dv) moves faster and experiences a lower energy barrier. The activation energy for
the first jump of Dv (0.06eV/b) is very similar to the one corresponding to a rigid
dislocation but after the maximum, the energy does not go down to the initial relaxed
value. As explained below, this is because of the presence of a kink pair.

Ta
double kink formation
22

c::

Q.,

;:,-.

150

200

z position (A)

(a)
1.6 .--.-----,----,---,--------,----,---,--------,--~~
II
II
II
II
II

t = 37.5 ps

1.3

100
z position (A)

(b)

150

200

Chapter4

Figure 4-10. (a) Profile of dislocation Dv (y position of the dislocation along the
dislocation line) at different times form 20 to 37.5 ps. (b) Core energy along the
dislocation line for the same times.

In Figure 4-lO(a), we show the position profile of dislocation Dv (z and y
represent the position of the dislocation in [111] and [11-2] directions in the plot) from 20
ps to 37 .5 ps in the MD simulation. Figure 4- lO(b) shows the core energy along the
dislocation line at the same times. We can see that at t = 20 ps the dislocation is almost
perfectly straight and the core energy is constant along the dislocation line. Note that the
position of the dislocation and core energy are affected by the presence of the vacancy
with a z position only close to zero. At time t = 30 ps a kink pair can be clearly seen both
from the dislocation profile and core energy plot (note that we have periodic boundary
conditions). Part of the dislocation line has advanced to the next equilibrium position
while the middle part of Dv (from z "" 60A to z "" 125A) is still climbing the Peierls
potential barrier. From Figure 4- lO(b) we see that the core energy is lower for the
portions of the dislocation that advanced and is still higher in the middle part of our
simulation cell. The asymmetry in the energy plots comes from the fact that the kink pair
contains two different kinds of kinks (will be explained in the next chapter).

4.7 Conclusion

Using the first principles based qEAM force field, we studied the mobility of the
a/2<111> screw dislocation in Ta. Applying the definition that the core of the l/2a
screw dislocation is formed by the 12 atoms with higher strain energy per Burgers vector,

Chapter4

100

we examined the variations of the core energy as the dislocation migrates. This leads to a
novel way of calculating the Peierls energy barrier and stress from MD simulations that
provides detailed information about the mobility of dislocations. This method gives a
Peierls energy barrier of Ep = 0.032 eV/b for twinning shear and Ep = 0.068 eV/b for
anti-twinning shear. The predicted Peierls stress at OK is ;=790 MPa for twinning shear
and ;=1430 MPa for anti-twinning shear. These values are about 7% larger than the
Critical Resolved Shear Stress (CRSS) calculated (using the same force field) by shearing
a large cylinder containing a dislocation. As in experiments and earlier simulations, we
find a clear non-Schmid behavior. The analysis of atomic strain energy distribution also
allows us to follow the path of dislocation migration under both twinning and antitwinning shear.

We find that both in the twinning and anti-twinning motion, the

dislocations move in <112> directions on {110} planes, but the actual path taken by the
dislocations differ due to the twinning/anti-twinning asymmetry of the energy landscape.

Our simulations at temperatures T = 20K, SOK, IO0K, and 300K have shown a
marked difference from those performed at T=0.00lK revealing the importance of
temperature effect on dislocation mobility. At high temperatures, thermal fluctuations
lead to incoherent motions of the segments within the dislocation apparently aiding the
migration. This leads to activation energy of 0.0053 eV/b.

Based on this activation

energy, we estimate that the Peierls flow stress for temperatures in the range 20K-300K is

~ 120 MPa. This is in good agreement with experimental results 1 1, 170 MPa at T=73K to
50 MPa at T =273K. Furthermore, in Section 4.6, our preliminary simulation on the kink
migration process in screw dislocations clearly shows that the screw dislocation may
decrease the energy barrier that impede its motion by forming a kink pair along its line. A

Chapter4

101

detailed study on the kink formation energy, structure and relations will be given m
Chapter 5.

4.8 Appendix: Periodic boundary for the simulation cell containing a screw dislocation
dipole
Our computations simulate the system as a crystal with periodic boundary
conditions, since this removes questions of the boundary surfaces and simplifies the
calculations. We construct the periodic cells containing a dipole of screw dislocations by
starting with a perfect periodic crystal and applying isotropic elastic theory. This causes a
partial stacking fault along the periodic boundary for the crystal cell unless the lattice
parameters for the crystal cell are optimized.

102

Chapter4

-2

-1

-------r---------1

-1:

r---------r---------r---------r---------r--------nL2

-2:

L---------L---------L---------L---------L--------nLt---..-...1
I•
nLl - - -......- L t • I •

Figure 4-11. The 2-D Schematic map of an atomistic simulation cell (solid primary
rectangle cell) containing a screw dislocation dipole and its n (n=2) layers of image cells.
The lattice parameters of the simulation cell are L 1 and L 2 . The Burger vectors of the
screw dislocation are normal to the plane. In the cells, the dislocation with Burgers vector
b (represented by a plus sign) is at the fractional coordinate (1/2, 3/4) but the dislocation

with Burgers vector -b (represented by a minus sign) is at the fractional coordinate (1/2,
1/4). The atom A and A' are on the boundaries parallel to the dislocation dipole and
equivalent in the periodic perfect crystal. While the atom B and B' are on the boundaries
perpendicular to the dislocation dipole and also equivalent in the perfect crystal.

Chapter4

103

Figure 4-11 shows the scheme we used to compute the atomic displacements for a
periodic screw dislocation dipole from isotropic elastic theory. The rectangular primary
dipole cell (in the center of Figure 4-11) is surrounded with n layers of its periodic image
cells. The Burgers vectors of the screw dislocations in the dipole are normal to the plane,
leading to atomic displacements only along the direction normal to the plane. Each cell
contains the dislocation dipole, where the dislocation with positive (pointing out) Burgers
vector has fractional coordinates of (1/2, 3/4) and the dislocation with negative (pointing
inside) Burgers vector is at (1/2, 1/4 ). The displacements for the atoms in the primary cell
are calculated by summing the contributions from all dislocations in the supercell, which
includes the primary cell and image cells. The calculated atomic displacements approach
their converged values as the number of image cells is increased.

We will consider next the issue of periodic boundaries m the calculations.
Consider two cases:

Equivalent atoms A and A' are on the boundaries parallel to the dislocation dipole
in the crystal cell before the introduction of the screw dislocation dipole.

Equivalent atoms B and B' are on the boundaries perpendicular to the dislocation
dipole.

We will now evaluate the displacement difference between these pairs of atoms
caused by the periodic screw dislocation dipole.

Figure 4-11 shows that the displacement of the atom A caused by the positive
dislocation in the image cell (1, 0) is same as the displacement of the atom A' caused by
the dislocation with the same Burgers vector in the image cell (1, 1). This is because the

Chapter4

104

lines (solid) connecting the atoms to the dislocations are parallel. However, there is no
dislocation in the supercell causes the same amount of displacement to the atom A as the
dislocation in the leftmost column of the image cells displaces the atom A' (as indicated
by the dashed lines). Thus the displacement difference between the atom A and A' is the
displacement of the atom A caused by the dislocations in the rightmost image cells less
the displacement of the atom A' caused by the dislocations in the leftmost column of the
image cells, as given in Eq. (4).

b n
(i + 3_ )L2 - X
(i + _!_ )L2 - X
13.dA-A' = -(L)tan- [
]- tan- [
]} +
Jr i=O
(n+-)L
(n+-)L1

(4)

(i+'}_)L2 +x
(i+_!_)L2 +x
}) tan -1 [
4 1
] - tan -1 [
4 1
] })
i=O
(n+-)L1
(n+-)L1
n-l

Here L1 and L2 are the lattice parameters of the primary cell, x is the distance from
the atom A (or A') to the bottom of the primary cell and n is the number of layers of the
image cells.

A similar procedure leads to the displacement between atoms B and B' given in
Eq. (5).

13.d B-B' = Jr

(i +_!_)LI - X

cI {tan
i=O

-l [

(n+_!_)L 2

n-1

L {tan
i=O

-1 [

(i + _!_ )L1 - X

] - tan -l [

]} +

(n+'}_)Lz

(i+_!_)Ll +x
(i+_!_)Ll +x
23
] - tan -1 [
2 1
]} )
(n + - )Lz
(n + - )L1

(5)

105

Chapter4

Here xis the distance from the atom B (or B') to the left boundary of the primary cell.

As n goes to infinity, we find that

f..dA-A' =-tan
1[

/'id 8 _ 8 , = 0,

-1

L2
[-],
LI

for all XE [0,L2 ]

(6)

for all x E [O, L,]

(7)

Eq. (6) implies that for finite values of L1 and L2 , introducing the periodic screw
dislocation dipole makes atoms A and A' nonequivalent. As a result, a partial stacking
fault along the boundary parallel to the dislocation dipole is formed in the crystal cell
with the magnitude of this stacking fault determined by the ratio of lattice parameters.
This stacking fault disappears only when L 2 << L 1, which corresponds to dislocation
dipole annihilation when L 1 is finite.
On the other hand, Eq. (7) shows that the screw dislocation dipole does not cause
a stacking fault along the boundary perpendicular to the dislocation dipole in the crystal
cell.

Chapter4

106

4.9 References
1. J.P. Hirth and J. Lothe, Theory of dislocations (Wiley, New York, 1982), p.241
2. G. Schoeck, Philos. Mag. A, 79, 2629 (1999).
3. L. H. Yang, P. Soderlind, and J. A. Moriarty, Philos. Mag. A, 81, 1355 (2001).
4. K. Ito and V. Vitek, Philos. Mag. A, 81, 1387 (2001).
5. M. S. Duesbery and V. Vitek, Acta Mater., 46, 1481 (1998).
6. M. S. Duesbery, Proc. R. Soc. Lond. A, 392, 175 (1984).
7. M. Parrinello and A. Rahman, J. Appl. Phys., 52, 7182 (1981).
8. W. G. Hoover, Phys. Rev. A, 31, 1695 (1985).
9. A. Strachan, T. (_;agm, 0. Giilseren, S. Mukherjee, R. E. Cohen, and W. A. Goddard,

First Principles Force Field for Metallic Tantalum, Phys. Rev. B, submitted.
10. D. E. Segall, T. A. Arias, A. Strachan, and W. A. Goddard, Accurate Calculations of

the Peierls Stress in Small Periodic Cells, J. Comput.-Aided Mater. Design, to be
published.
11. R. Lachenmann and H. Schultz, Scripta Met., 4, 709 (1970).
12. V. Vitek, Cryst. Lattice Defects, 5, 1 (1974).
13. G. Wang, A. Strachan, T. Cagin and W. A. Goddard, Journal of Material Engineering
and Science A, 309-310, 133 (2001).

Chapters

107

Chapter 5 Flips and kinks on 1/2a screw dislocation in Ta

5.1 Overview
The plasticity of metals and semiconductors is controlled by the properties of
dislocations and the interactions between dislocations with other defects in crystals.
Hence, knowledge of the structure, self-energy, and evolution pattern of dislocations is
essential to arrive at a good understanding of plastic deformation of materials and to
obtain a mesoscopic model of deformation processes 1-4 • Much information on
dislocations can be obtained from such high-resolution experimental techniques as
HRTEM and STM. However, many details of the structural and energetic properties of
dislocations are beyond the resolution of current experimental methods. Computer
simulations at the atomistic level provide the best way to attain deeper insight about
dislocations 5 ·6 .

In bee metals (e.g., K, a-Fe, Mo, and Ta) at low temperatures, the crystal lattice

resists the motion of screw dislocations more strongly than the motion of edge
dislocations 7 . Therefore, the mobility of screw dislocations governs the plastic
deformation behavior of these materials at low temperatures. From atomistic simulations
at OK, the screw dislocation is thought to move in a rigid, collective fashion leading to a
minimal external Peierls stress of about 10- 2 µ (µ is the shear modulus of the crystal) 8- 11 •
However, the observed rapid decrease of the Peierls stress with increasing temperature
implies that at finite temperatures the screw dislocations move by formation and
subsequent migration of kinks rather than by translation of the straight dislocation 12 .

108

Chapters

The concept of kinks in dislocations and using kinks in describing plastic flow
behavior of crystal were mathematically treated in the framework of elasticity theory by
Seeger and Schiller 13 in 1966. These ideas are still applicable. The first direct observation
of the dislocation kinks was made by Kolar et al. 14 using atomic resolution transmission
electron microscopy (TEM) on partial dislocations in Si. Many modern mesoscale
plasticity theories (for instance, Ref. 4) use the kink-pair mechanism to describe the
motions of dislocations. These theoretical models can benefit from the accurate atomistic
descriptions of dislocation kinks provided in this thesis.
Using an atomistic simulation, Seeger et al. 15 proposed that the asymmetric
dislocation cores for the 1/2a screw dislocation in a-Fe were polarized, and then
explained the multiplicity of kinks and the existence of flips (antiphase defect5) in
dislocations. In two classical papers 16 •17 , Duesbery studied the detailed structure, Peierls
stress, and formation energy of the isolated kinks in the 1/2a screw dislocation in K
and a-Fe. Duesbery and Basinski 18 showed that atomistic computer simulations of kink
pair generation and migration agreed with the experimental flow stress of Potassium (K).
Recently, the formation energies of kinks in screw dislocation in Ta 9 and Mo 19 have
been determined much more accurately in simulations with Green's function boundary.

In this chapter, we use a simulation model with periodic/ fixed boundaries to
(1) determine the formation energies of 1/3a<112> kinks in the 1/2a screw

dislocation in Ta,
(2) estimate the lateral motion energy barriers of those kinks,
(3) analyze the configurations and structures of thosekinks,

Chapter5

109

(4) investigate the inherent relationship between different types of kinks.

The remainder of this chapter is organized as follows. Section 5.2 describes of the
details the periodic/fixed boundary simulation models. Section 5.3 describes the core
configurations of the 1/2a screw dislocation and outline the types of the dislocation
defect (flip and kink). Section 5.4 reports our results on formation energies of the isolated
kinks and kink pairs, while Section 5.5 estimates the migration energies of kinks without
using an applied stress. Section 5.6 describes our analysis of the flip and kink structures
and the inherent relationship between different kinks. In this section, we also summarize
and explain the trend of the kink formation energy and mobility from present work and
literature 19 . Finally, our conclusions are given in Section 5.7.

5.2 Simulation model

5.2.1 Construction of simulation model

To study kinks in dislocations, we use the model crystal shown schematically in
Figure 5-1, which is orthorhombic and oriented by the [11-2], [l-10] and [111] crystal
directions. This cell consists of three distinct construction regions (region A, region B
and region C) in the [111] direction.

The construction region A and region C contain four a/2<111> screw dislocations
arranged as a quadrupole, in which a pair of dislocations has Burgers vector b = a/2[111]
and the other pair of dislocations has Burgers vector b = a/2[-1-1-1]. In region A and C,
the initial dislocation was constructed based on elasticity theory and it was subsequently

110

Chapters

relaxed to reach its equilibrated configuration under 3-D periodic boundary conditions.
The positions of the dislocations in the region A and region C differ by a vector v from
the equilibrium dislocation center in the region A pointing to the equilibrium dislocation
center in the region C as indicated in Figure 5-1.

----------------/

---- __________ y ------------/

A------i~B-+1<11----- C

U1211f [11 O]

4r 111 1
Figure 5-1.

The schematic plot of the simulation model. In this model, region A and

region C contain the equilibrated dislocation quadruples. Region B is constructed based
on elastic theory to smooth the configuration misfit. The vector v starts from the center
of the dislocation in region A and points to the dislocation in the region C. In our
simulations, v can only be O (flips), 1/3a[ll-2] (right kinks) or 1/3a[-l-12] (left kinks).
The shaded regions indicate the fixed boundaries, which are 5 b thick, in the simulation.
The cell parameters are 5a[ll-2] (=40.7 A), 9a[l-10] (=42.3A), and 150a/2[111]
C=43L8A).

111

Chapters

The central region B is designed to smooth the interfacial misfit between the
region A and region C. The initial atomic displacement relative to the perfect crystal for
each atom in the region B is obtained from elastic theory by the following equations.
(1)
with
r·[lll]

a=--hB

(2)

where ~d[;111 (r) and ~d[~111 (r) are the displacements determined by elastic theory for the
atom positioned at r caused by the periodic dislocation quadruples in the region A and
the region C, respectively. The h8 is the height of the region B in the [111] direction.

5.2.2. Boundary conditions of simulation model
We employed the qEAM many-body force field (FF) 20 to describe the atomic
interaction potentials for Ta. This embedded-atom-model force field was derived from
accurate quantum mechanics (QM) calculations. It describes with good accuracy the bee,
fee and A15 phases of Ta for pressure from ~ -10 0Pa to ~ 500 0Pa, and also the
vacancy formation energy, surface energy and shear twinning energy for bee Ta crystal.
The qEAM FF has previously used to study the melting temperature of Ta as a function
of pressure 20 , spall failure 21 and properties of straight dislocations 22 .

The simulations impose periodic boundary conditions in the [11-2] and [1-10]
directions of the model crystal. The quadruple arrangement of the 1/2a screw
dislocation in the (111) plane eliminates the misfit of atoms on the periodic boundaries
due to dislocation images. On both ends of the simulation cell along the [111] direction,

Chapters

112

there are fixed regions that are 5 b (~ 14.4 A) [larger than the cut-off radius (9 A) of the
qEAM FF] thick. In this way, a 3-D simulation model is formed without introducing
misfit or free surface. The movable atoms interacting with the fixed boundaries in the
simulation effectively interact with an infinite equilibrium dislocation quadrupole and do
not "feel" the existence of free surface.

It is appropriate for us to discuss the effects of the boundary conditions and the
size of the simulation cell on the resultant kink at this point.

1) First, the 2-D periodic boundary of the model makes the simulation easy to
implement. However, it introduces arrays of the dislocations with the kinks in the
(111) plane. Also there is interaction energy between the kinks on different
dislocations. This kink-kink interaction energy must be taken into account when we
calculate the formation energy of the isolated kinks. Section 5.4 uses isotropic
elasticity theory to estimate the kink-kink interaction energy raised by the 2-D
periodic boundary and corrects to obtain the formation energy of the isolated kink.

2) Second, the fixed boundary in the [111] direction might cause an atomic misfit near
the boundary if proper caution is not exercised. We computed the final kink width
(details see Section 5.6.4) and kink formation energy of a NRP kink (definition see
Section 5.3) using simulation cells with different lengths in the [111] direction. The
results (from line 1 to line 4) in Table 5-1 demonstrate that the kink formation energy
as well as the kink structure is well converged when the simulation cells are more
than 100 b long. The equilibrium kink width obtained from the strain energy
distribution converges to 18 b, while the geometric kink width converges to 10.7 b. In

113

Chapters

a simulation cell of length 38 b, incomplete relaxation causes the calculated kink with
a 22% wider strain energy peak and 12% larger kink formation energy.
3) The third issue is the height of the region B in Figure 5-1. The optimal choice is the
width of the kink of interest. However this information is not available before the
simulation. The analysis in Section 5.6 shows that 10 b is a good approximation for
all kind of kinks. In fact, our results (line 4, 5 and 6 in Table 5-1) show that the choice
of this height has no affect on the kink structure and the kink formation energy.
Table 5-1. The determined final kink width and kink formation energy (before the
correction of the kink-kink interaction) for a NRP kink (definition see Section 5.3) in
different simulation cells. All simulation cells are 40.7 A long in the [11-2] direction and
42.3 A long in the [1-10] direction. In the table, the total length of the simulation cell and
the length of the region A, B and C as indicated in Figure 5-1 are given in the unit of
Burgers vector b, which is 2.88 A.

Length in the [111] direction (b)

Final kink width (b)

Formation
energy (eV)

Total

Region A

Region B

Region C

Method I

Method II

9.9

221

0.696

10.7

0.625

136

10.7

0.624

150

10.7

0.624

150

10.8

0.624

150

10.7

0.624

Chapters

114

With the above considerations, we employed the simulation cell whose geometry
was 5a[l-12] (=40.7 A), 9a[l-10] (=42.3 A) and 150a/2[111] (=431.8 A) in our study. As
indicated in Figure 5-1, the length of the region A, region B and region C is 70 b (=201.6
A), 10 b (=28.8 A) and 70 b (=201.6 A), respectively. Our cell contains 40,500 atoms

(37,800 movable) in the simulation cell.

5.3 Multiplicity of flips and kinks

5.3.1 Equilibrium dislocation core structure

A. Differential displacement map

We used elasticity theory to construct the initial simulation cell with screw
dislocation quadrupole. Then, we used the qEAM FF to minimize the total strain energy
of the quadruple, obtaining the equilibrium dislocation configuration.
The differential displacement (DD) map 23 in Figure 5-2 (left column) shows the
local strain field around the dislocation center. The interpretation of these DD maps is
given in the figure caption. Figures 5-2 (a) and (b) show two equilibrium dislocation
cores for the dislocation with b is equal to 1/2a[ 111] while Figures 5-2 (c) and (d) show
the dislocation cores when b is 1/2a[-l-1-1]. These figures show that the equilibrium
dislocation core has threefold symmetry and spreads out in three <112> directions on the
{ 110} planes in the DD map. There are 6 equivalent <112> directions on the (111) plane,
so there exist two kinds of core configurations both for the dislocation with Burgers
vector 1/2a[ll l] and for the dislocation with Burgers vector 1/2a[-1-1-1]. Despite the

Chapters

115

difference in core configurations, all four dislocations have the same energy both in terms
of self-energy (including core energy and elastic energy) and in terms of elastic
interaction energy. Note that the quadrupole has four dislocations and each of which has
one of two possible core configurations. This causes 16 combinations of the quadrupole,
but we find that all have identical energies.

B. Relaxation map

The displacement of each atom along the [111] direction from the atomistic
relaxation as compared to that calculated from continuum elasticity theory is shown in
the relaxation maps (the right column of Figure 5-2). The magnitude of this difference for
all atoms, except the 6 columns of atoms closest to the dislocation line, is less than 0.05

A (0.017 b). These atoms are at distance more than 1.44 b (4.15 A.) from the center of
dislocation. This demonstrates that elasticity theory describes the elastic field of screw
dislocation quite well and fails only near the core region of the dislocation (within 1.09 b
= 3.31 A.). The direction and magnitude of the displacement difference for the central 6

columns of atoms are presented in the relaxation maps. The most important result in these
maps is that three central atoms of the dislocation relax simultaneously 0.267 A (=0.09 b)
either in the [111] direction in a P (Positive) type dislocation or in the [-1-1-1] direction
in an N (Negative) type dislocation. This phenomenon is called the polarization of
dislocation 15 . Regardless of the orientation of Burgers vector b, the P type dislocation
core spreads along the [-1 -1 2], [-1 2 -1], and [2 -1 -1] directions in the DD map. While,
the N type dislocation cores spread out along the [1 1 -2], [1 -2 1], and [-2 -1 -1]
directions.

Chapters

(a) N+

116

• •

o-•-•
/\

• •

-o

(b) P+

• •

• •
•-•--o
I\
' oI
\ I
• •
,-·,
• •
• •
• I \

...:1:-·

l...,j

o-•-•
• •
• •
,,
• •
I \ •

(d) P-

'.)

. n:
o-

I \

•-•-o
• •
• •

• •

-0.267

o.m

o.Lti.o •
-0.267

-0.267

0.123

• •

• 0Tl"'" l T0•

-O.m •

-O.m

0.267
0.267
-0.123

• •

• •
• •
• • t • •

()

0.123
-0.267
-0.267

,., r• 0•
• 0•r. T
;·"'
o,.m.

• •

-0.123

• 0·1r1·0 •

"·
• T
T •
-0.123

-0.123

0.267

• •

Chapters

Figure 5-2.

117

The equilibrated dislocation core configurations for the 1/2a screw

dislocation in Ta. The circles represent the projected atoms in the (111) plane. The open,
shaded or black circles indicate that the atoms are in three consecutive ( 111) layers of bee
lattice. However, the arrows in two columns of figures have different meanings.

The left column shows the differential displacement map, in which the arrow
indicates the displacement in [111] direction (perpendicular to the map) of the
neighboring atoms relative to their positions in the perfect bee crystal. The direction of
the arrow represents the sign of the displacement and the magnitude is proportional to the
relative displacement between corresponding atoms. When the arrow touches the two
atoms, the relative displacement between these two atoms is 1/3 b. For clarity, the
relative displacements less than 1/12 bare not shown in the figure.

The right column shows the relaxation map, in which the arrow from each atom
indicates the relaxation (parallel to the dislocation line) relative to the displacement field

predicted by isotropic elastic theory. The magnitudes of such relaxation (in angstrom) for
the central 6 columns of atoms (the relaxation for the other atoms is less than 0.05A) are
printed next to the corresponding atom. Four types of energy degenerate dislocation core
configurations are distinguished in terms of the relaxation direction of the three central
columns of atoms (downward denoted as "N" and upward denoted as "P") and Burgers
vector (a/2[111] denoted as "+" and -a/2[111] denoted as "-").

Chapters

118

5.3.2 Flips
By definition, the flip 24 (or antiphase defect) is a defect where a core
configuration changes to the other one along the screw dislocation line. Two kinds of the
1/2a screw dislocation (N-type and P-type) lead to two possible configurations of
flips (from P to N and from N to P) as schematically shown in Figure 5-3(a). We find the
formation energy of a P-N flip is 0.005 eV suggesting this flip could occur thermally
along the dislocation. The nucleation energy of the N-P flip is 0.572 eV, which is higher
than the energy fluctuations. The P-N and N-P are two distinct flip configurations in
l/2a screw dislocation.

5.3.3 Isolated kinks

The kink refers to the region in which one segment of the dislocation in an energy
minimum connects with another segment that lies in a neighboring position. In this study,
we focused our interest on the kinks where the dislocation segments are separated by
either 1/3a[ll-2] (called the Right kinks) or -l/3a[ll-2] (called the Left kinks). Figure 53(b) shows that in each category (Right or Left) of the kinks there are four combinations
of the dislocation core configurations. We thus have 8 possible kinks: NRP, NRN, PRP,
PRN, NLP, NLN, PLP and PLN. Among these kinks, the NRN and PRP are energy
degenerate and related by symmetry operations, so are the NLN and PLP.

Chapter 5

(a)

119

Flips

P-N
Llli=0.005 e V

(b)

N-P
Llli=0.572 eV

Right kinks ( v =_!_a[l 12])

NRP

NRN

PRP

PRN

8E=0.655 eV

Llli=0.634 eV

8E=0.634eV

8E=0.611 eV

Left kinks ( v = _!.a[TT 21 )

NLP

NLN

PLP

PLN

Llli=l.153 eV

8E=0.632 eV

Llli=0.632 eV

8E=0.139 eV

Figure 5-3. The schematic drawing, nomenclature and calculated formation energy of the
defect (flip and kink) in 1/2a screw dislocation. In figures, the arrow (.._) represents
P type dislocation and ~) represents N type dislocation. (a) Two kinds of flips exist in

screw dislocation. The core configuration along a straight dislocation line can flip either

Chapters

120

from P to N (denoted as P-N) or from N to P (denoted as N-P). The formation energy of
N-P (0.572 eV) is larger than that of P-N (0.005 eV). (b) There are 4 kinds of right kinks
(NRP, NRN, PRP, and PRN) and 4 kinds of left kinks (NLP, NLN, PLP, and PLN). The
defect vector v (indicated in Figure 5-1) is 1/3a[ll-2] for right kinks and 1/3a[-1-12] for
left kinks.

5.4 Kink formation energy calculations

5.4.1 Formation energy of the isolated kink

With the knowledge of the dislocation core configurations, the dislocation
polarization and the multiplicity of flip and kink in the 1/2a screw dislocation, it is
feasible to construct a dislocation quadruple such that all dislocations in it have the same
type of defects (flip or kink). We constructed and relaxed the quadrupole of the
dislocations with the defect in the way described in Section 5.2. After obtaining the
equilibrated dislocations with the defect, we calculated the total energy of the relaxed cell
[Ed(cell)] by summing the atomistic energy for all movable atoms in the simulation. This

energy includes the self-energies of the dislocations with the defect [Ed( self)] and the
interaction energy between the dislocations with the defect [Ed(inter)].

On the other hand, the total energy [Ep(cell)] of a quadrupole with the same
geometry but containing four equilibrium straight 1/2a screw dislocations can also
be calculated with simple 3-D periodic boundary simulation (the fixed boundaries in
Figure 5-1 are removed). Similarly, the Ep(cell) can be expressed as the self-energies of

121

Chapters

the dislocations without the defect [Ep( self)] plus the interaction energy between the
perfect dislocations [Ep(inter)].
The formation energy of the defect (flip or kink) is the self-energy difference
between an isolated dislocation with the defect and the dislocation without that defect

Thus the intrinsic formation energy of a defect (iJ.E1) is expressed as

M 1 = : [Ed (cell) - E P(cell)] - : [Ed (inter)- E P(inter)],

(3)

In Eq. (3), the first term _!_[Ed(cell)-EP(cell)] (called the differential cell energy) is

obtained directly from the simulations while the second term _ _!_[Ed(inter)-EP(inter)]

(called the interaction correction) is obtained from elasticity theory.

In elastic theory, the flip in the dislocation is considered as a dimensionless point
defect and has no effect on the interaction energy between dislocations. Thus the second
term is Eq. (3) is O when calculating the flip formation energy. However, this term does
not vanish in the case of kink, which is a dislocation defect with finite dimension. The
interaction energy between two kinked dislocations and the interaction energy between
two straight dislocations can be calculated by summing the contributions from all
piecewise straight segments 26 . This approach has been used to derive the elastic energy of
the kink pair in the same dislocation 25 . Furthermore, the converged value of the second
term in Eq. (3) can be obtained by summing the pair interactions in the 2-D periodic
quadrupole of the kinked dislocations and the straight dislocations.

We consider two simple models for the description of kinks.

Chapters

122

The perpendicular model assumes the kink is a pure edge segment, which is
2.71 A (ll/3a1) long in the <112> direction. For every kind of kink in
the l/2a screw dislocation, isotropic elastic theory comes to 0.030 eV
for the second term of Eq. (3).

The inclined model. In fact, the equilibrated kink is not a line segment
perpendicular to the dislocation line but a region whose height is about 2. 71 A
in the <112> direction and whose width is around 10 b (~28.8 A) along the
[ 111] direction (see Figure 5-14 ). The inclined model assumes the kink is a
dislocation line segment spanning a width w along the dislocation line and a
height h normal to the dislocation (the values of wand h of the kink are given
in Table 5-5 of Section 5.6.4). Assuming isotropic Ta, The energy difference
between a pair of kinked dislocations and a pair of straight dislocations,
denoted as W(L 1, L 2 ), is calculated using the following equations.

(4c)

123

Chapters

(4d)

µb 1b 2
h L~ ·[/(Li,L2 )+Rw(Li,Li)]+
h w L; ·I(Li,L2 )
]'
4Jr(l-v)[
h L~ +w (L~ +L;)
(w +h )[h L~ +w (L: +L;)]

In the above equations, L1 and L2 are the separation distances between
dislocations in the [11-2] and [1-10] directions; wand hare the kink width and height; b 1
and b2 are the Burgers vectors of two dislocations. The shear modulus µ is equal to C44
and the Poisson ratio v = 2C12/(C 11 +C12 ). In the calculations, we adopted the elastic
moduli (C11=272.54 0Pa, C12=137.57 0Pa and C44=69.63 0Pa) of the perfect bee Ta
crystal at OK determined by the qEAM FF20 .

Table 5-2 gives the calculated results for the differential cell energy [first term of
Eq. (3)] from the simulation and the interaction correction [the second term of Eq. (3)]
from the inclined model calculation, which we used to calculate the formation energies
for various flips and kinks. The values of the interaction correction from the inclined
model deviate by at most 0.003 eV from the 0.030 e V obtained assuming the kink
perpendicular model. This implies that ignoring the real geometry of the kink causes only
a marginal error in determining the formation energy (e.g., 0.1 % for the NRP kink and 2
% for the PLN kink). Thus the calculated kink formation energy is insensitive to the

uncertainty of the kink geometry.

Table 5-2. The differential cell energies (e V) from the qEAM FF simulations, interaction
corrections (eV) from continuum theory using the inclined model [Eq. (4)], and the

124

Chapters

intrinsic formation energies (eV) of the defects (flips and single isolated kinks) in the
1/2a screw dislocation in Ta.

Differential cell energy a

Interaction correction b

Intrinsic formation energy

I/4[Eicell)-Ep(cell)]

-I/4[Eiinter)-Ep(inter)]

AEr

N-P (flip)

0.572

P-N (flip)

0.005

NRP (right kink)

0.624

0.031

0.655

NRN (right kink)

0.604

0.030

0.634

PRP (right kink)

0.604

0.030

0.634

PRN (right kink)

0.582

0.029

0.611

NLP (left kink)

1.122

0.031

1.153

NLN (left kink)

0.601

0.031

0.632

PLP (left kink)

0.601

0.031

0.632

PLN (left kink)

0.106

0.033

0.139

configuration

a The perpendicular model gives 0.030 eV.

b see Eq. (3).

5.4.2 Formation energy of kink pairs

In addition to a single kink at which the dislocation line crosses a Peierls energy
hill, there are also kink pairs consisting of a left kink and a right kink. These kink pairs
can be formed by thermal fluctuation in the crystal and their nucleation and subsequent
motion are thought to be important in low temperature deformation processes of bee
metals. If the separation between the left and right kink is sufficiently large, the formation

Chapters

125

energy of a kink pair is just the sum of the formation energies of the two component
kinks. Since there are 4 kinds of left kinks and 4 kinds of right kinks, there are 16 ways to
combine a pair of the kinks in the 1/2a screw dislocation in Ta. In some cases, one
or two flips are required to fulfill the requirement of the dislocation core configuration
when the kink pair nucleates from a perfect dislocation.
The whole spectrum of the configurations and the formation energies of all
possible kink pairs are given in Table 5-3. The calculated formation energies of the kink
pairs range from 0.794eV to 1.894 eV. We find that the PLN-NRP kink pair has the
lowest formation energy, which is 0.794 eV. This formation energy is close to the value
of 0.81 eV for the zero shear stress activation enthalpy of the 1/2a screw
dislocation in Ta determined by Tang et al. 27 by fitting the empirical data to the Kocks
model. Our calculated range covers the available experimental measurements (0.92 eV
by Funk28 ,1.24 eV by Rodrian et al. 29 , 0.98 eV by Wemer30 , and 0.97 eV by Mizubayashi
et al. 31 ) of the formation enthalpy of the double-kink on the 1/2a screw dislocation

in the Ta single crystal. Our calculated kink pair formation energy compare favorable
with of 0.88eV to 1.50 eV calculated by Yang et al. 9 using the multi-ion interatomic
potential from the model generalized pseudopotential theory (MGPT) for Ta. This
agreement is somewhat surprising. In our study, the three columns of the atoms closest to
the dislocation core shift 0.09 b along the dislocation line causing an asymmetric
dislocation core. In the MGPT FF calculations these atoms only translate 0.0007 bin the
[111] direction leading to a symmetric core. The agreement between these two
calculations suggests that the symmetry of the l/2a screw dislocation core does not
play an important role in the study of screw dislocation kinks.

126

Chapters

Table 5-3. Calculated formation energies of all kink pairs in 1/2a screw dislocation
in Ta. The total formation energy of kink pair is the summation of the formation energies
of the component single kinks and the flips required. Note that the kink pair PLN-NRP
has the lowest formation energy, which is 0.472 eV lower than the second lowest kink
pair formation energy.

configuration

Initial flip

Right kink

Left kink

Formation energy (e V)

NRP

PLN

0.794

NRN

NLN

1.266

PRP

PLP

1.266

NLN

1.292

NRP

Internal flip

P-N

NRP

PLP

1.292

N-P

PRP

PLN

1.345

PLN

1.345

NLP

1.764

NLP

1.792

NRN

N-P

PRN
PRP

P-N

NRN

NLP

1.792

N-P

PRN

NLN

1.815

PRN

N-P

PLP

1.815

P-N

NRP

P-N

NLP

1.818

N-P

PRP

P-N

NLN

1.834

P-N

NRN

N-P

PLP

1.834

N-P

PRN

N-P

PLN

1.894

Chapters

127

5.5 Kink migration energy calculations

5.5.1 Kink migration energy
Once the kink pair nucleates, the component kinks would move laterally driven
by an applied resolved shear stress. During the lateral motion, the kink would experience
periodic energy barriers from crystal lattice. If the required kink migration energy were
comparable with the kink formation energy, both kink pair formation and migration
processes would govern the mobility of the dislocation. Hence, it is also important to
quantify the kink migration energy. In this section, we propose a way to estimate the
magnitude of the kink migration energy and demonstrate the difference of motion for
various kinds of kinks.
In the simulation cell containing the equilibrated dislocations, the position of each
atom differs from its position in perfect bee crystal by an amount of !ir, which is the
atomistic displacement. If there is no kink in the dislocation, the atoms in the same
column in the [111] direction will have the exactly same atomistic displacement.
However, the existence of the kink in the dislocation destroys such regularity. The atoms
in the kink region have different atomistic displacements from those atoms far away from
the kink region even though they are in the same [11 l] column. When the kink migrates
along the dislocation line one step, the strain field of the whole simulation model will
migrate along the [ 111] direction by 1b, as will the atomistic displacements. In the
current study, we translate the strain field rigidly and estimate the energy barrier during
the kink moves along the dislocation. Suppose two consecutively neighboring atoms,
atom i and atom j, are in a [111] column and the corresponding atomistic displacements

Chapters

are

128

!1r 0 and !1rJ0 in an equilibrium dislocation with a kink. After the kink moves a

distanced in the direction from the atom i to the atom j, the atomistic displacement of the
atom i is determined with the following equation.
A -h
A -o
A -o
LJ.r.
= (1 - -d)
. LJ.r.
+ -d . LJ.r.
J '

(5)

A new configuration representing the moving kink is obtained by updating the atomistic
displacements for all atoms in the model crystal.
We calculated the potential energy for every configuration and determined the
potential energy barrier as the kink moves one Burgers vector along the dislocation line.
The Sr/ in Eq. (5) keeps unchanged for a perfect dislocation because/1,; 0 is equal to

11rf , so the perfect dislocation segments do not cause any variation in the potential
energy. Therefore, the calculated energy barrier must be the kink lateral migration
energy. The kink migration energy of a NRP kink is estimated to be 2.5 x 10-4 eV (0.04%
of its formation energy 0.655 eV) and the PLN kink migration energy is 3.5 x 10-4 eV
(0.3% of its formation energy 0.139 eV). The NLP kink was found to have the largest
migration energy 1.9 x 10-3 eV, which is only 0.2% of its formation energy 1.153 eV. In
our calculation the kink moves in a rigid and collective way, which makes our result an
overestimate of the kink migration energy. It applies only in the limit of low stress
deformation. Since the kink migration energy is about two magnitudes smaller than the
corresponding kink formation energy, it is evident that at low stress conditions the
mobility of the 1/2a screw dislocation in Ta is controlled by the kink pair
formation energy rather than kink migration energy. The same conclusion has also been
drawn from MGPT FF simulations 9 • Because the kink migration energy is much less

Chapters

129

significant than the kink formation energy, we do not think an accurate determination of
the kink migration energy is a priority.

5.5.2

Relative mobility of kinks

In section 5.3 and section 5.4, we studied the multiplicity of kinks and computed
the formation energies of different kinks. It is also interesting to investigate the kink
migration mobility for different kink configurations. The method described in subsection
5.5.1 can be used to estimate the kink migration energy barrier well. However, this
method does not emphasize the role of the kink configuration in the migration process
because atoms in the kink region and in regions far from the kink are translated
simultaneously. In the following, we propose a way to compare the migration mobility of
kinks.

Analyzing the atomic motions during kink migration, we found that atoms in the
kink region move much more than atoms far from this region. Hence, we can partition the
atoms in the simulation cell into two groups (group A and group B) using a cutoff
parameter y. When kink migrates by 1 b along the dislocation, the atoms in group A
translate more than y while the atoms in group B move less than or equal toy. Using this
grouping strategy, we constructed the dislocation configurations describing the kink
migration process as follows: the atoms in group B are positioned as they were in the
equilibrium kink configuration while the positions of the atoms in group A are computed
by Eq. (5). In fact, we believe that this trajectory is close to what occurs to the kink when
it moves rapidly under a high stress.

Chapters

130

We carried out one-point potential energy evaluations for the different
configurations of kink motion. The equilibrium configuration of the kink has the lowest
potential energy and the potential energy of the system increases during the kink
migration. The energy increase, called the elastic relaxation energy, indicates the far field
of a kink (composed of the atoms in the group B) resists the migration of the core region
of kink. Using a cutoff parameter y as 0.05 A, Figure 5-5 shows the elastic relaxation
energies for different kink configurations [all right kinks in 5-5(a) and all left kinks in 55(b )]. The internal friction between the atoms in the group A and those in the group B
causes the elastic relaxation energy increasing quadratically. The elastic relaxation
energy after the kink moves 1 b can be used to infer the mobility of the kink (assuming
that a higher elastic relaxation energy implies it is much harder to move the core region
of the kink along the dislocation line).

131

Chapters

0.14

Right kinks in Ta

G-E) NRP kink

1 = 0.05t\

G--EJ NRN kink
~ PRPkink

l:r--t:, PRN kink

0.02

Reaction coordinate

(a)

0.14

Left kink~ in Ta
G-E) PLN kink

y"' 0.051\.

G--EJ PLP kiuk

'1-----v NLN kink
l:r--t:, NLP kink

.~
er,

~ 0.04
{J.02

Rcartion coordiuatc

(b)

Chapters

132

Figure 5-5. The elastic relaxation energy associated with the kink lateral migration by +lb. Atoms that move less than or equal to 0.05A between the initial and final equilibrium
configurations are kept fixed as their initial positions. While the atoms, which move more
than 0.05 A between two equilibrated configurations, are moved rigidly by 0.1 b at each
step by linear extrapolation. In our computation, there are 14 such atoms per NRP kink,
30 atoms per NRN or PRP kink, 43 atoms per PRN kink, 15 atoms per PLN kink, 25
atoms per PLP or NLN kink and 36 atoms per NLP kink. (a) Right kinks and (b) Left
kinks. The results show that the NRP kink and the PLN kink have the lowest elastic
relaxation energy when kink move 1 b such that they have the highest migration mobility
among the right and the left kinks, respectively.

For the various kinds of kinks, we plot in Figure 5-6 [all right kinks in Figure 56(a) and all left kinks in Figure 5-6(b)] the elastic relaxation energy after the kink moves
1 b. In these figures, we show the calculated elastic relaxation energy under various cutoff parameter y. The results indicate that the mobility of kinks differ appreciatively when

y > 0.01 A but are close when y < 0.01 A. Therefore, we expect to observe a mobility
difference between kinks in high-stress conditions but similar migration behavior for the
kinks under low stress. For y = 0.05 A, the calculated the elastic relaxation energies when
kinks moved 1 bare about 0.1 eV. This hints that the kink migration energy would play
an important part under high stress conditions. These computations show that the
migration mobility of kinks is in the order of:
For right kinks, NRP > NRN (=PRP) > PRN,

(6a)

For left kinks, PLN > NLN (=PLN) > NLP,

(6b)

An explanation of the relations in (6) will be given in the next section.

133

Chapters

0.14

Right kinks in Ta

G-f) NRP kink
[3-£] NRN kink

9-s.:/ PRP kink
i:s:--6 PRN 1.i nk

-~
] 0.08

...

,tj

:::
::)

.§ 0.04

:E"'
0.02

{Hll

0.02

0.03

0.04

0.05

(a)

0.14

Lel't kinks in Ta

G-f) PLN kink
[3-£] PLP kink

9-'v NLN kink
i:s:--6 NLP kink

1 0.06

"i:i
:::

.§ 0.04

i"'

0.02

0.03

(b)

0.04

(J.05

134

Chapters

Figure 5-6. The maximum elastic relaxation energy (when kink moved 1 b) associated to
the kink lateral migration. The calculations same as reported in Figure 5-5 have been
carried out for the different y (from 0.001 to 0.05 A). (a) Right kinks and (b) Left kinks.

The above results allow us to conclude that among all possible kink pairs the
PLN-NRP kink pair has not only the lowest formation energy but also the lowest
migration energy barrier.

5.6 Structural analysis

5.6.1 Overview

In this section, we present a detailed structural analysis of dislocation defects (flip
and kink). The study aims to elucidate the reasons for the following:
(1) Why do N-P and P-N flip have different formation energies?
(2) Why does formation energy of the NRP, NRN (or PRP) and PRN kinks
decrease and differ by ~0.02 eV while the formation energies of the NLP,
NLN (or PLP) and PLN kinks decrease but differ by ~0.50 eV?
(3) Why does the mobility of kinks follow the rule: NRP > NRN (=PRP) > PRN
and PLN > NLN (=PLN) > NLP?

We also carried out structural analyses on kinks to determine the geometrical
parameters of the isolated kink (kink width w and kink height h) and to estimate the
minimum stable separation between a pair of kinks.

Chapters

135

5.6.2 Structural analysis of flips
A. Relative displacement between neighboring columns (DD map)

There are two kinds of flips for the a/2<111> screw dislocation in Ta. They are
the P-N flip, whose formation energy is 0.005 eV, and the N-P flip, whose formation
energy is 0.572 eV. Figure 5-7(a) and Figure 5-7(b) show the strain energy distribution of
the relaxed quadrupole of dislocations containing the P-N flips or the N-P flips along the
dislocation line. The strain energy is computed by summing the atomistic strain energies
(the atomistic energy for each atom in the simulation cell less the atomistic cohesive
energy in the perfect bee Ta crystal) for all atoms in a lb thick slice region. For the sake
of comparison, the strain energy distribution of a perfect dislocation quadrupole in the
same size simulation cell is also plotted. Figure 5-7(a) and Figure 5-7(b) show that the
strain energy of the dislocations with flips deviates from that of the perfect dislocations
only in the flip formation region (30 b long for P-N and 50 b long for N-P). It is
interesting that the middle 10 b (Z from 70 b to 80 b) of a dislocation with P-N flip
[Figure 5-7(a)] has less strain energy than the perfect dislocation while one with N-P flip
in Figure 5-7(b) has a strain energy maximum. The DD maps for both flips in Figure 57(c) show the atomistic configurations of the dislocation core at various positions marked
in Figure 5-7(a) and Figure 5-7(b). The DD maps for the P-N and N-P flip at the center of
flip [figure B and E of Figure 5-7(c)] are extremely similar. At the flip center, the
dislocation core is symmetric with zero polarization, quite different from the equilibrated
dislocation cores [shown in the figure A, C, D, and F of Figure 5-7(c)].

136

Chapters

7.9
.'iE = 0.005 eV/llip

~ with P-N flip

7.8

ca--o without l1ip
A: l'•IJp< dnlocation iZ: 40 b >

7.7

B, flip from P to N iZ: 75 b.i
C: N•typ<• di'llorntion ii.= 110 b)

1.)

A [Fig. .5-7(cl]

7.4

7.3
B [Fig. .5-7(c)]
7.2

100

110

Z fb)

(a)
7.9
AE = 0.572 cV/flip

~ with N-P flip

7.8

->
.......

..c

'.l
;;,.,

...

,:,J)

7.7

ca--o

without flip

D: N•IYJJ<' dislorn!lm1 iZ =40 b)

£, flipfromN foP!/.:=75 b)

E [Fig. 5-7(c.1]

F: l'•l~p• dislocation I Z=llO b)

7.6

1.)

,;::

·c::::....

7.5

7.4

{;,

E-

7.3

7.2
40

80
Z (b)

(b)

100

Chapters

137

P-N

N-P

ri'. ..1.,...

e .- G.,. 0

• •

B~:. • .. •

• '

• •

• ,•\ ZS:I•' ,•

• \ •. .\ 0....

.......
/\

....... 0

• '

• •

(c)

• •

138

Chapters

Figure 5-7.

The strain energy distribution for the dislocation quadrupole with (a) P-N

flip (formation energy is 0.005 eV) and (b) N-P flip (formation energy is 0.572 eV). (c)
The differential displacement maps show different core configurations along the
dislocation line (Z=40 b, 75 b and 100 b). Note: Although the dislocation with the P-N
flip has a lower strain energy than the perfect dislocation in the P-N flip formation region,
the total strain energy of the dislocation with a P-N flip is still 0.005 eV higher than the
total strain energy of perfect dislocation.

B. Relative displacement within a column

The DD maps in Figure 5-7(c) show the relative displacements in the [111]
direction between the neighboring atoms in a (111) plane in the dislocation core region.
They do not contain information on the relative displacements in the [111] direction
between the neighboring atoms in the same [111] column. In a perfect bee crystal or a
crystal with a straight 1/2a screw dislocation, the distance between two
neighboring atoms in the same [111] column is 1 b (ll/2a1). However, because of
the change of the polarization along the dislocation, this regular atomic separation is
expected to change when a flip is formed in the 1/2a screw dislocation. For each
individual column of atoms in the dislocation core, we calculated the distances between
two consecutively neighboring atoms in the [111] direction. Figure 5-8 shows the atomic
arrangement in a dislocation core. In this figure, the circle in the (111) plane represents a
[111] column of atoms and the cross mark indicates the center of the dislocation. The
atoms in the columns marked with the same letter are energetically and geometrically
equivalent by symmetry.

139

Chapters

00000

0@®®©0
©®®®@
0®®®®0
0®®®0
00@©00
Figure 5-8. The (111) projection of atom arrangements around an a/2[111] screw
dislocation. The cross indicates the center of the dislocation. The atoms marked with the
same letters are energetically equivalent and related by symmetry. For the P type
dislocation core configuration, the atomistic strain energies are in the order of the
A>B>C>E>D>F; while for the N type dislocation core configuration, the order is
A>B>D>F>C>E.

Figure 5-9(a) shows for the P-N flip the distance in the [111] direction between
consecutive atoms in the same [111] column varies along the dislocation and Figure 59(b) shows the same plot for the N-P flip. In both cases, the distance between
neighboring atoms in the [111] direction is equal to 1 b when far from the flip formation
region but deviates significantly within the flip formation region. Figure 5-9(a) shows
that for the P-N flip the distance between "A" atoms is compressed to 0.976 b (2.81 A)

Chapters

140

while Figure 5-9(b) shows that for N-P flip it is stretched to 1.032 b (2.97 A). In
comparison, for the P-N flip the "B" atoms are in tension (the maximal distance is
1.010b) while for the N-P flip they are in compression (the minimal distance is 0.996 b).
The atoms "C", "D", "E", and "F" also have different mechanical states but with smaller
magnitudes.

C. Energy distribution along the dislocation line

Figure 5-9(c) shows for the P-N flip the atomistic strain energy for each atom
along the dislocation while Figure 5-9(d) shows the same for the N-P flip. Obviously, the
change of the distance between neighboring atoms in the [ 111] direction affects the
atomistic energy for atom in the flip formation region. Summarizing the atomistic strain
energies for the 18 marked atoms comes to the core energy [shown in the insets of Figure
5-9(c) and 5-9(d)]. The elastic energy for each flip is computed by subtracting the core
strain energy contribution from the total strain energy of a dislocation containing a flip.
Figures 5-9(e) and 5-9(f) show the calculated elastic energy distribution along the
dislocation with P-N or N-P flip.

The middle 10 b at the center of a P-N flip has higher elastic strain than the
perfect dislocation energy but lower core energy. In comparison, the N-P flip has lower
elastic energy in the middle 10 b but a higher core energy than the perfect dislocation.
The energy field of a flip is the region where the energy (elastic or core energy) deviates
from that of a perfect dislocation. The elastic energy field for a flip of 30 b for P-N flip
and 50 b for N-P flip is much longer than its core energy field (20 b for P-N flip and 15 b
for N-P flip). For a P-N flip (formation energy is 0.005 eV), the core strain energy part is

Chapters

141

-0.275 eV while the elastic energy part is 0.280 eV. Similarly, the formation energy of
the N-P flip (0.572 eV) can be partitioned in to 0.773 eV core energy and -0.201 eV
elastic energy.

1.03

1.02

P-N flip

G-€J atom A

atom B
t:r--i:,. atom C and D
~ atom E and F
[3---£]

= L0l
!>Li

'?.
"';.)
(.)

1.010 b !BJ

1------

1.002 b [E. F]

~ 0.9'>

0.998 b [ C, DJ

0.98

0.97
20b
0.96 O

80
Z (b)

(a)

100

110

142

Chapters

L04
L03

1.032 b (A]

G---€J atom A
G--EJ at om 13

N-P flip

tr-i:l atom C and D

t.02

"?--'v atomE and F

,£.
:=' LOI

1.003 b [C. D]

§l•
1:3

:::
,,. <199

0.996 b LE, F]

0.98

0.980 b [BJ
1-e

0.97
0.96

~ I

15 b

100

110

Z (b)

(b)

0.3

1.65

P-N flip

core strain energy (e\/·

1.6
1.55

0.25

1.45

·\o 50 60 70 80 90 100 0

0.2

atomA

..::,

;.;.....

,...at,

g 0.15
1.J

ato111 B

___

(/l

0.1

,.... atomD

atomC

/~----

atomF

0.05

70
Z tb)

(c)

100

110

143

Chapters

1.65 core strain energy (cV
0.3

N-P flip

1.6
l.55

1.5

0.25

3 0.2

atomA

,,,,:,...,,

,,g 0.15

aton1 B

(/)

0.1

aton1 D

0.05

atomC

~o~F==,
50

atomE

- - 80

100

Z (b)

(d)
0.38
@-€) with P-N flip

- - pcrfcrt dislocatiou

0.37

.g 0.36
0::

"Vi

ti 0.35

0.33

30b

0. 32 ---~~-~~--'-----'----"'---'----'---'------'--'-----'-----'
40
50
60
70
80
90
100
llO
Z (b)

(e)

Chapters

144

0. 36 r----r----.-.-----,----.-,----,------,~-,------,-r----.----,------,
CH) with N-P flip
- - perfect dislocation

0.35

.§ 0.34

:-a

ti 0.33
c:...

,1,;

-~ 0.32

0.31

50b
0.3

100

1To

Z (b)
(t)

Figure 5-9. (a) The distance between neighboring atoms in the same [111] column around
an a/2<111> screw dislocation with a P-N flip. (b) The same plot as (a) for an N-P flip.
(c) The strain energy by atom around the 1/2a screw dislocation with a P-N flip.
(d) The same plot as (c) for an N-P flip. The positions of the atoms relative to the
dislocation center are shown in Figure 5-8 using the same letters. The insets in (c) and (d)
plot the summation of atomistic strain energy for individual atoms. Note that the obtained
core energy for the equilibrium dislocation is different from calculations in Ref. 22,
where the core energy is defined as the summation of 12 atoms with highest atomistic
strain energy. To study the strain energy change in the flip region, we include 18 atoms in
this computation. (e) The elastic energy [one-fourth of total energy in (a) less core strain
energy in (c)] along a dislocation with P-N flip. (f) The same plot for a dislocation with
N-P flip.

Chapters

145

D. Discussion
The different ways to flip the polarization of dislocation cause that the same atom
is at different mechanical states (compression or tension) in the [111] direction and has
the different atomistic strain energy. This is the reason why the formation energy of P-N
flip is different from that of N-P flip.

The calculated formation energy of N-P flip is 0.23 eV in Ref. 9, 0.20 eV in Ref.
19 and 0.572 eV in the present work. While the formation energy of P-N flip is 0.03 eV
in Ref. 9, 0.00 eV in Ref. 19 and 0.005 eV in the present work. Our analysis shows that
the atoms B (see Figure 5-8) in the different mechanical states are the principle cause that
the formation energy of the N-P flip is higher that that of the P-N flip. The atoms B are
stretched along the [ 111] direction and contribute -0.183 eV to the P-N flip formation
energy. However, these atoms are compressed in the [111] direction and give 1.242 eV to
the N-P flip formation energy in our study.

It is also accountable that the formation energies of two kinds of flips in

1/2a screw dislocation will differ with a less magnitude (0.20 eV in Ref. 9
compared to 0.567 eV in our study) when the dislocation cores are only slightly
polarized. The smaller polarization difference of different types of dislocations implies
less compression or tension in the [111] direction for atoms.

146

Chapters

5.6.3 Structural analysis of kinks

A. Relation of kinks

The formation energies of kinks given in Figure 5-3 clearly show the following
trends. For right kinks, the formation energies decrease in the order of NRP, NRN
(=PRP), and PRN with differences of ~0.02 eV. For left kinks, the formation energies
increase in the order of PLN, NLN (=PLP), and NLP with differences of ~0.5 eV. We
carried the structural analysis of these kinks to understand the origin of these trends.

Figure 5-10 shows the strain energy distribution maps for various right kinks.
These figures show that
(1) The NRP kink [Figure 5-lO(a)] has only a strain energy maximum at its

formation region,
(2) the NRN kink [Figure 5-lO(b)] has a strain energy maximum at the formation
region and a strain energy minimum above its formation region,
(3) the PRP kink [Figure 5-lO(c)] has a strain energy maximum at the kink
formation region and a strain energy minimum below the formation region,
(4) the PRN kink [Figure 5-lO(d)] has a strain energy maximum at the kink
formation region and two strain energy minima on both sides of the formation
region.

147

Chapters

7.9
AE = 0.655 eV/kink

G-€) with NRP kink

E (Fig. 5-1 ll

7.8

13--!J

without kink

A: Equilibrium S t~pf' disloratiou

E: Tramition "ifat(> of 1'RP I.ink.
ll: f.uilibrium P 1ype dislorntion

A(f'ig. 5-Ul
·.-~

n (Fig. 5-UI

7.4

r7.3

7.'2

100

110

Z (b)

(a)
7.9
AE = 0.634 eV/kink

G-€) with NRN kink

7.8

.- 7.7

LJ
_...,
►,

13--!J

without kink

A: Euilibrium -1' 1y~ cllidocation

E: Tr311Stion tilate of NRP link
F: Trmnition ,late of P•l\ flip

7.6

O.L

....

C: f.rntion

,t;
C:
,:I

...
'iii

7.5
A (Fig. 5-11 I

7.4

t-'

7.3

7.'2

Z (b)

(b)

100

110

Chapters

148

7.9

ll.E = 0.634 cV/kink

G-0 with PRP kink
G-EJ without kink

7.8

E (Fig. 5-111
7.7

>u
~·
>-. 7.6

'-<

D: Equilibrium P l~p< dislocation
C: Traruiition ,tat< of P•l\ Oip
E~ l'.-amilion s1ate of SRP h.ink
Jl: E.qiulibrium P typ< dislocation

7.5

D (Fig. 5-111

8 (Fig. 5-1 l 1

::;

7.4

t""'

7.3
G (Fig. 5-Hl

7.2
40

100

1fo

Z fb)

ll.E = 0.611 cV/kink
7.8

G-0 with PRN kink
c.--o

E (Fig. 5-HI

without kink

D: EcJUilibrimn P tnie dnlocation

7.7

G: Transition '!it.ate of P--N Dip
E, Transition slate of :-.RP I.ink

F: Tramrition .tale of I'•N flip
C: :Equilibrium N type dislorntfon

D (Fig. 5-11.1
c::

7.4
7.3

F (Fig. 5-111

7.2
40

70
Z (b)

(d)

100

110

Chapters

Figure 5-10.

149

The strain energy distribution for dislocation quadruples with right kinks.

(a) the kink, (b) the NRN kink, (c) the PRP kink, and (d) the PRN kink. The letters in the
figures indicate the regions with characteristic features along the dislocation. The
dislocation core configurations of these regions are shown in Figure 5-11.

Figure 5-11 shows the DD maps corresponding to these critical states. The panels
(A), (B), (C), and (D) in Figure 5-11 show the dislocation core configurations at regions
far from the kink formation region. These maps are same as those of the equilibrated
dislocation cores. The configurations of dislocation core in the region where the strain
energy is a maximum in all four right kinks have the same differential displacement
pattern as shown in Figure 5-ll(E), which clearly indicates that the whole dislocation
evenly splits into two parts in the neighboring equilibrium positions. We find that Figure
5-1 l(F) and (G) resemble the DD map of the flip in the panel C of Figure 5-5(c) such that
the strain energy minima in Figure 5-10 correspond to the flips in screw dislocation.

150

Chapters

~···········································~

•'

~----------------------------------------····
.... o

c::;

• •

....

• •
• •

• •

...... 0

,,
I\'
.. .

• •

.......••....•......................••..•••

r•••••••••••••••••••••••••••••••••••••••••••~

p:::;

O•O

• 0 •

r•••••••••••••••••••••••••••••••••••••••••••~

• • '
~-------------------------------------------~

,•

• ..I

~-··········································~
-~. . . . . . o •• :

••O

:w•,.

\ I

O••••

e • 0 e :
.••••••••••••••••............................

Chapters

151

Figure 5-11. The differential displacement maps of dislocation core at different regions
marked in Figure 5-10 along the l/2a[ll 1] screw dislocation. The figures (A), (B), (C),
and (D) show the equilibrium dislocation cores; the figure (E) shows the atomic relative
displacements at the center of the kink formation region; while the figures (F) and (G)
indicate the flips in the kink formation region.

On the above analysis, the relation of the right kinks in the 1/2a screw
dislocation can be expressed as the following equations.
NRN = NRP + P-N,

(7a)

PRP = P-N + NRP,

(7b)

PRN = P-N + NRP + P-N,

(7c)

These equations indicate that the NRP kink is the elementary right kink and all other right
kinks are composites consisting of the NRP kink and one or two P-N flips. The NRP kink
and the P-N flips are only separated by 3 bin the composite kinks. The formation energy
of an isolated P-N flip is 0.005 eV, but the close distance between the NRP kink and the
P-N flip may decrease the total strain energy and leads to a -0.02 e V contribution for each
P-N flip in the composite kink. Thus, Eq. (7) also explains why the formation energy of
the NRP kink is 0.021 eV higher than that of the NRN (or PRP) kink and 0.044 eV higher
than that of the PRN kink.

A similar analysis for the left kinks is presented in Figure 5-12 and Figure 5-13.
Figure 5-12 shows the strain energy distributions along the dislocation with the different
left kink.

152

Chapters

7.9
AE = 0.139 eV/kink

Go-€) with PLN kink
r.---a

7.8
7.7

without kink

A: £

£: Trmwition 'Yt:ite of PL"'li lJ:111'.
ll: Equih1>rium :0.. IJJ"' di,locatinn

B (Fig. 5-131

A (Fig. 5-131

7.3
E (Fig. 5-131

7.2
40

100

110

Z (b)

(a)
7.9
AE = 0.632 cV/kink

7.8

;,e
>~
;:,-,

7.7
A: £

7.6

,../;Jj
·;;

F: 'I rnnliffion ,tale of :,...p flip
C: F.

7.5

...

F (Fig. 5-13.1

E: Trmwition slate of PL"/ kink

A (Fig. 5-13)

C (Fig. 5-131

7.4

f-,

7.3
7.2
40

80
Z rb)

(b)

100

110

153

Chapters

7.9
AE = 0.<>32 eV/kink

@-€) with NLN kink

7.8

[r-1]

I>: Equilibrium ~ type disloc-ation

7.7

G: Transition ~fate ofN--P flip

"-'

...

without kink

E: Tnm,rition .tale of PL", IJ.nb

G (Fig. 5-131

7.6

H: Equilibrium '.\' type di'Uncation

-~1.

-; 7.5

-;;;

D (Fig. 5-BI

B ( Fig. 5- Bl

7.4
7.3

7.'2

110

100

z tb)
(c)
7.9
7.8

AE = 1.153 eV/kink
E (Fig. 5-13)
D: F..quilib rium N type dislocation
C: Trnnfiilfon ,tale ofN•P Olp
E, J'nuwi!ion .tale of PLN kin!.

F: Tnm,ilion •late of:\'•P Oip
C: Equilihdum p type ruslocatioll

F (Fig. 5-BI
G (Fig. 5-Bl
D ( Fig. 5-13)

C (Fig. 5-131

7.3

7.'2
40

80
Z (h)

(d)

100

110

Chapters

154

Figure 5-12. The strain energy distribution for the dislocation quadruples with left kinks.
(a) the PLN kink, (b) the PLP kink, (c) the NLN kink, and (d) the NLP kink. The letters
in the figures indicate the regions with characteristic features along the dislocation. The
dislocation core configurations of these regions are shown in Figure 5-13.

In these figures, we see a strain energy minimum at the PLN kink formation
region and a superficial resemblance of the strain energy distribution, in which there is
only a strain energy maximum at the kink formation region, for the NLN, PLP and NLP
kink. It seems that there is no obvious relationship between the left kinks. However, we
are still able to establish the linkage between left kinks scrutinizing the detailed
differential displacement maps in Figure 5-13. Figures 5-12(F) and (G) strongly suggest
the existence of the N-P flip in the formation region of the NLN, PLP and NLP kink.
Thus, the PLN kink is the basic left kink. All other left kinks are the combinations of the
PLN kink and one or two N-P flips as indicated in the following equations.
NLN = N-P + PLN,

(Sa)

PLP = PLN + N-P,

(Sb)

NLP = N-P + PLN + N-P,

(Sc)

The formation energy of an isolated N-P flip is 0.572 eV. The above equations explain
the difference of 0.50 eV in the formation energies of the PLN (0.139 eV), PLP (or NLN,
0.632 eV) and NLP (1.153 eV) kink.

155

Chapters

.••......•.•••...........••.............•....
B(/ • • .. • 0

. . . . f;i;.•

• -

~Ml.
,~-------------------------------------····•·J

..... o

c~~

• I

'o\ ...
• •

• •

~----------------------------·········------J
fE:·; • •• 0 •

• •

• • I• \.. 0 •

...'.'

•\

•••

a:; • ,... o

,ow'\
\' • '
o ",
O+e
19

~--■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■--

,--------------------------------------------.,

... .....1, •
.... . •

• I\
• '•
,~-------------------------------------------J

,r•••••••••••••••••••••••••••••••••••••••••••.,
~J

\ I

.............................................

Chapters

156

Figure 5-13. The differential displacement maps of dislocation core at different regions
marked in Figure 5-12 along the 1/2a[l ll] screw dislocation. The figures (A), (B), (C),
and (D) show the equilibrium dislocation cores; the figure (E) shows the atomic relative
displacements at the center of the kink formation region; while the figures (F) and (G)
indicate the flips in the kink formation region.

Besides the relation of the kink formation energy, the relation of the kink mobility
in Eq. (6) can also be accounted by Eq. (7) and Eq. (8). The existence of flips in the kink
will decrease its mobility. Thus the NRP kink and the PLN kink have the highest mobility
in the right and left kinks.

B. Discussion

The kink relationship in Eq. (7) and Eq. (8) provides the first such connection in
the atomistic level simulations. Although these relations were obtained using a qEAM FF
for Ta, they provide a universal pattern for all bee metals. To prove this point, we
compared all available data of the kink formation energies in bee metals 9 •17 •19·32 .

A direct corollary of Eq. (7) is that the kink formation energy differences &NRN_
&NRP' &PRN_&NRN and (&PRN_&NRP)/2 should be nearly equal and close to &P-N'

which is the formation energy of the isolated P-N flip. Based on Eq. (8), the kink
formation energy differences &NLN_&PLN' &NLP_&NLN and (&NLP_&PLN)/2 should be
similarly close to the formation energy of the N-P flip (&N-P). It should be addressed that
the flip in the composite kinks (NRN, PRP, PRN, NLN, PLN, and NLP) is under the
different environments from the isolated flip. The close interaction between the flip and

157

Chapters

the kink might relax the total strain energy, such that the kink formation energy
differences could be smaller than the corresponding flip formation energy.
Table 5-3. Comparison of the formation energies (in eV) of the flips under different
environments.

a-Fe a

Mob

Tac

(Duesbery)

(Rao et al.)

(Yang et al.)

(present work)

~P-N

0.048

0.300

0.00

0.03

0.005

~NRN_~NRP

0.043

0.267

-0.16

-0.11

-0.021

~EPRN -~ENRN

-0.022

-0.085

-0.15

0.20

-0.023

l/2 (~PRN_~NRP)

0.011

0.091

-0.16

0.05

-0.022

~EN-P

0.018

0.408

0.21

0.23

0.572

~ENLN_~PLN

0.028

-0.322

0.18

0.19

0.493

~NLP_~NLN

0.045

0.126

0.21

0.08

0.521

l/2 (~ENLP _~PLN)

0.037

-0.098

0.20

0.14

0.507

Materials

P-N flip

N-P flip

a Reference 17, using a first-principle interatomic potential for potassium and an
empirical interatomic potential for iron.
b Reference 19, using the MGPT FF.
c Reference 9, using the MGPT FF.

Chapters

158

Table 5-3 compares the formation energy of the isolated flips and the flips in the
composite kinks. Both present results for Ta, and calculations by Rao et al. 19 for Mo
show the similar regularity of the flip formation energies as discussed above.
Furthermore, the plot of differential energy between the kinked dislocation and an
unkinked configuration for the PRP kink (denoted as p-pf kink in Figure 6 of Ref. 19)
resembles the Figure 5-lO(c) showing the strain energy distribution for the PRP kink in
present work. We believe that the energy minima in both figures indicate the existence of
a P-N flip in the PRP kink formation region. We used the qEAM FF for Ta as well as the
periodic boundaries in the [11-2] and [1-10] direction and the fixed boundary in the [111]
direction in our study. While the MGPT FF for Mo and Green's function boundary
conditions were employed in Ref. 19. The agreement between these two simulations
indicates that the relation of kinks in Eq. (7) and Eq. (8) is independent of the employed
force field and boundary conditions.
However, the results by Yang et al. 9 using the MGPT FF for Ta and Green's
function boundary conditions (also in Table 5-3) did not show a similar regularity of the
flip formation energies. Neither did the even older calculations by Duesbery 17 for K and
a-Fe. There are two possible reasons for this discrepancy. First, the equilibrium
dislocation core in our study has a large polarization (~ 0.09 b) whereas the dislocation
polarization is small (~ 0.0007 b) in Ref. 9. A smaller polarization of the dislocation
implies a smaller difference among the kinks in the same category (Left or Right). The
composite kinks might not dissociate into a flip and an elementary kink to decrease the
strain energy when the dislocation core is symmetric and only weakly polarized. The
second reason could be the incomplete relaxation of the atomistic structures. Duesbery

Chapters

159

used the fixed boundaries where atoms are fixed at the positions determined by
anisotropic elasticity theory in the simulation. These fixed boundaries could introduce
bias in the atomistic relaxation if the simulation cells were not sufficiently large in three
dimensions.

Ref. 32 found the following order of kink pair formation energies
PLN-NRP < NLN-NRN < NLP-PRN,

(9)

However, no atomistic explanation was proposed. It is easy for us to interpret Eq. (9)
with the help of the relation of kinks. The kink pair NLN-NRN can be considered as the
combinations of the kink pair PLN-NRP and a pair of the N-P and the P-N flips in the
kink region. Similarly, the kink pair NLP-PRN can be considered as the kink pair NLNNRN plus a pair of the N-P and P-N flips in the kinks. Such that if the pair of a N-P and
P-N flips contribute a positive strain energy in the composite kinks, the increasing order

in Eq. (9) would hold true. Actually, Eq. (9) is universal as demonstrated by Table 5-4.
All the available kink pair formation energies, except for a-Fe in Ref.17, follow the same
trend. As to the failure case, the empirical potential for iron yielded negative formation
energies for two kinds of kinks discredited those results. Thus, the kink pairs formation
energies in K, Mo, Ta, and a-Fe obey the rule (9), so far.

Table 5-4. Comparison of formation energies of kink pairs. In the table, "Yes/No"
indicates whether the calculated kink pair formation energies do or do not obey the rule:
PLN-NRP < NLN-NRN < NLP-PRN as in Eq. (9).

160

Chapters

~E PLN-NRP (eV)

~ NLN-NRN ( eV)

~E NLP-PRN (eV)

Yes/No

Ka (Duesbery)

0.076

0.147

0.170

Yes

a-Fe a (Duesbery)

0.241

0.186

0.227

Mo b (Rao et al.)

1.62

1.64

1.70

Yes

Tac (Yang et al.)

0.96

1.04

1.32

Yes

Ta (present work)

0.794

1.266

1.764

Yes

a-Fe ct (Wen et al.)

0.84

1.29

1.94

Yes

Materials

a Reference 17, using a first-principle interatomic potential for potassium and an

empirical interatomic potential for iron.
b Reference 19, using the MGPT FF.

cReference 9, using the MGPT FF.
ct Reference 32, using a nudged elastic band method and an EAM potential.

5.6.4 Determination of geometrical parameters

In addition to the formation energy and the migration energy of the kinks, the
geometrical parameters, such as the kink height hand kink width w, are also essential for
a mesoscopic description of kink in the continuum model. In this subsection, we present
our efforts to determine these parameters from the dislocation line shape and strain
energy distribution. The calculated results are given in Table 5-5.

Table 5-5. The computed geometrical parameters for the kinks in the 1/2<111> screw
dislocation in Ta. In Method I, the width of the kink is determined by fitting a straight
line to the kink formation region and determining the distance between the intersections
of this line with the two limited locations, which are two neighboring straight dislocation

161

Chapters

centers. In Method II, the width of the kink is determined to be the length of the region at
whose boundary the strain energy deviates by 0.2% from the equilibrium dislocation
strain energy.

Configuration

Height (A)

Width (b) [Method I]

Width (b) [Method II]

NRP (right kink)

2.77

10.7

NRN (right kink)

2.71

10.4

PRP (right kink)

2.71

10.4

PRN (right kink)

2.65

10.2

NLP (left kink)

2.65

8.9

NLN (left kink)

2.71

9.1

PLP (left kink)

2.71

9.1

PLN (left kink)

2.77

9.3

Figure 5-14 shows a line representing the dislocation with the (a) NRP, (b) NRN,
(c) PRP, (d) PRN, (e) PLN, (f) PLP, (g) NLN, and (h) NLP kink. Every point in the line
is determined by calculating the atomistic strain energy weighted center for those 12
atoms with the highest strain energy in a 1b slice of the dislocation. The figures show that
the dislocation is in its equilibrium position in the regions far away from the kink
formation region. The average distance between two equilibrium positions on two sides
of kink is the kink height. As indicated in Figure 5-14, the heights of the kinks are not
equal for different combinations of dislocation core configurations.

Chapters

162

2.5
,.-.,

i:-<

2.5

(H).'.'iRI'l.lnk

('~ 1.5

2.771'.

(b)
G-E) :-R'.'i kink

Nl.5
......

0.5

40 50 60 70 80 90 100 I 10
Z (b)

40 50 60 72 (b~O 90 100 110

2.5

(d)

<,:::,
N 1.5

N l.5

0.5

40 50 60 70 80 90 HXJ lfo
Z (b)

40 50 60 70 80 90 100 110

(e)

2.5

-- 2

N'l.5

0.5

0.5
50 60 70 80 90 100 110

50 60 70 80 90 HXJ 110

Z (b)

(g)

2.5

r,,<

(f)

2.5

o< 2

Z (b)

--<<

Nl.5

(h)

2.5

NL5

- 1

2.(i5J\

0.5

0.5
50 60 70 80 90 100 110
Z (b)

40 50 60 70 80 90 HXJ 110
Z (b)

Chapters

163

Figure 5-14. The profile of the dislocation lines in the kink formation regions. (a) NRP,
(b) NRN, (c) PRP, (d) PRN, (e) PLN, (f) PLP, (g) NLN, and (h) NLP kink.

The kink width w can be estimated in two ways. (I) The part of the dislocation
line in the kink formation region (70b :S Z :S 80b shown in Figure 5-14) was fitted into a
straight line. The kink width is the distance in the [111] direction (the abscissa in Figure
5-14) between two intersections of this line with two equilibrium dislocation lines
separated by the kink height. (II) The kink width is the length of the region bounded by
two points, where the strain energy deviates from the strain energy of the equilibrium
dislocation by 0.2%, in the dislocation.

Determining the kink width by the line shape of the dislocation in Method I and
the strain energy distribution in Method II is suitable for different applications. The kink
width from Method I is the geometrical description of the kink and was used to compute
the kink-kink interaction energy in Section 5.3. In the mesoscale model 4 , the minimum
stable distance between a left kink and a right kink is required. Instead of carrying out
simulations checking the stability of the kink pair positioned at various separations, we
can estimate the minimum stable distance between kinks using the kink width determined
in Method II. The strain energy distributions of a pair of kinks do not overlap each other
at their minimum stable distance, so the minimum stable separation of this pair of kinks is
one half of the summation of two component kink widths determined in Method II. The
choice of 0.2% in the calculation is somehow arbitrary. However, the attained minimum
stable distance of the NRP-PLN kink pair is 16 b, which is close to the 13 b obtained by
an empirical fit 4 . So, 0.2 % is a reasonable choice.

Chapter 5

164

5.7 Conclusion
In this paper, we report calculations on the kink formation energies and the
equilibrium kink structures using the first principle qEAM FF. The formation energies of
the kink pairs in 1/2a screw dislocation are in the range of 0.794-1.894 eV, which
agrees with the results of the MGPT FF calculations and is in the same range of the
empirical data. The PLN-NRP kink pair was found to have the lowest formation energy
that is 0.794 eV compared favorably with the 0.81 eV for the zero shear stress activation
enthalpy from the empirical data fitting. Our detailed structural analysis reveals that the
PLN kink and the NRP kink are the elementary left and right kinks. The other kinks are
the composite kinks composed of the elementary kink and the flips. This relation of kinks
accounts for the observed trend of the kink formation energies and mobility.

As an atomistic simulation, our results are limited by the description capability of
the force field. Our simulations yield the asymmetric and polarized (~0.09 b) dislocation
cores, which cause the multiplicity of kinks and the existence of flips. However, the
recent ab initio calculations 33 ' 34 obtain only a symmetric screw dislocation core in Ta,
also simulations using the Finnis-Sinclair type atomic interaction potential come to a
symmetric core 35 . The MGPT FF calculations 9 yield an almost symmetric (slightly
polarized ~0.0007 b) dislocation core in Ta. It is clear that dislocation polarization is very
sensitive to the calculations. Some dislocation properties depend on whether the
dislocation core is symmetric or asymmetric. For instance, the formation energy
difference between the N-P flip and the P-N flip in this study (dislocation polarization is
0.09 b) is 0.567 eV. It is more than two times larger than the 0.20 eV, when the
dislocation core is symmetric and polarized only by 0.0007 b. However, some other

Chapters

165

dislocation properties would not be very sensitive to the difference of the dislocation
cores. We showed in a preceding paper

that the discrepancy of the obtained dislocation

cores did not affect the calculated core energy and Peierls stress. In this paper, our results
of kink formation energies are consistent with the results from the MGPT FF, though the
symmetry of dislocation core is different in two studies.

On the other hand, the

polarization of screw dislocation is subject to a rapid change depending on the volume
and pressure conditions as pointed out by Yang et al. 9 The physics, such as the relation of
kinks proposed in present work, should be still applicable when the dislocation is
polarized under certain conditions. So, the difference of the obtained dislocation core
configuration between the present work and other studies should not undermine the
credibility of the present work.

Chapters

166

5.8 References

1. V. V. Bulatov and L. P. Kubin, Curr. Opin. Solid State Mater. Sci., 3, 558 (1998).
2. R. Phillips, D. Rodney, V. Shenoy, E. Tadmor, and M. Ortiz, Model. Simul. Mater.
Sci. Eng., 7, 769 (1999).
3. M. I. Baskes, Curr. Opin. Solid State Mater. Sci., 4, 273 (1999).
4. A. M. Cuitino, L. Stainier, G. Wang, A. Strachan, T. <;agm, W. A. Goddard, and M.
Ortiz, A Multiscale Approach for Modeling Crystalline Solids, J. Comput.-Aided
Mater. Design, to be published.
5. V. V. Bulatov, S. Yip and A. S. Argon, Philos. Mag. A, 72,453 (1995).
6. J. A. Moriarty, W. Xu, P. Soderlind, J. Belak, L. H. Yang and J. Zhu, J. Eng. Mater.
Tech., 121, 120 (1999).
7. F. Louchet and L. P. Kubin, Philos. Mag. A, 39,433 (1979).
8. W. Xu and J. A. Moriarty, Comput. Mat. Sci., 9, 348 (1998).
9. L. H. Yang, P. Soderlind, and J. A. Moriarty, Philos. Mag. A, 81, 1355 (2001).
10. K. Ito and V. Vitek, Philos. Mag. A, 81, 1387 (2001).
11. G. Wang, A. Strachan, T. <;agm, and W. A. Goddard, Mater. Sci. and Eng. A, 309,
133 (2001).
12. T. Suzuki, Y. Kamimura and H. 0. K. Kirchner, Philos. Mag. A, 79, 1629 (1999).
13. A. Seeger and P. Schiller, Physical Acoustics, edited by W. P. Mason (Academics,
New York 1966), Vol. 3A, p. 361.
14. H. R. Kolar, J.C. H. Spence and H. Alexander, Phys. Rev. Lett., 77, 4031 (1996).
15. A. Seeger and C. Wuthrich, Nuovo Cimento, 33B, 38 (1976).

Chapters

167

16. M. S. Duesbery, Acta Metall., 31, 1747 (1983).
17. M. S. Duesbery, Acta Metall., 31, 1759 (1983).
18. M. S. Duesbery and Z. S. Basinski, Acta Metall. Mater., 41, 643 (1993).
19. S. I. Rao and C. Woodward, Philos. Mag. A, 81, 1317 (2001).
20. A. Strachan, T. <;agm, 0. Gtilseren, S. Mukherjee, R. E. Cohen, and W. A. Goddard,

First Principles Force Field for Metallic Tantalum, Phys. Rev. B, submitted.
21. A. Strachan, T. <;agm and W. A. Goddard, Phys. Rev. B, 63, 060103 (2001).
22. G. Wang, A. Strachan, T. <;agm, and W. A. Goddard, Atomistic Characterization of

the Screw Dislocation in Tantalum, Phys. Rev. B, submitted.
23. V. Vitek, Cryst. Lattice Defects, 5, 1 (1974).
24. J. P. Hirth, Acta Mater., 48, 93 (2000).
25. J. P. Hirth and J. Lathe, Theory of dislocations (Krieger, Melbourne, FL, 1982),
p.242.
26. J. P. Hirth and J. Lathe, Theory of dislocations (Krieger, Melbourne, FL, 1982),
p.162.
27. M. Tang, L. P. Kubin and G. R. Canova, Acta Mater., 46, 3221 (1998).
28. G. Funk, Dissertation, University Stuttgart, 1985.
29. U. Rodrian and H. Schultz, Z. Metallk, 73, 21 (1982).
30. M. Werner, Phys. Stat. Sol. (a), 104, 63 (1987).
31. H. Mizubayashi, H. Egashira and S. Okuda, Acta Metall. Mater., 43,269 (1995).
32. W. Wen and A.H. W. Ngan, Acta Mater., 48, 4255 (2000).
33. S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett., 84, 1499 (2000).
34. C. Woodward and S. I. Rao, Philos. Mag. A, 81, 1305 (2001).

Chapters

168

35. M. S. Duesbery and V. Vitek, Acta Mater., 46, 1481 (1998).

Chapter 6

169

Chapter 6 A multiscale approach for modeling crystalline solids *
6.1 Overview
The proposed multiscale modeling approach for advanced materials (such as highpurity bee single crystals) is aligned with the current divide and conquer paradigm in
micromechanics 1-6 . This paradigm first identifies and models the controlling unit
processes at microscopic scale. Then, the energetics and dynamics of these mechanisms
are quantified by means of atomistic modeling. Finally, the macroscopic driving force is
correlated to macroscopic response via microscopic modeling. This last step involves two
stages, localization of the macroscopic driving force into unit-process driving forces and

averaging of the contribution of each unit process into the macroscopic response.
We will show that the meticulous application of this paradigm renders truly
predictive models of the mechanical behavior of complex systems. In particular, we
predict the hardening of Ta single crystal and its dependency for a wide range of
temperatures and strain rates. The feat of this approach is that predictions from these
atomistically informed models recover most of the macroscopic characteristic features of
the available experimental data, without a priori knowledge of such experimental tests.
This approach then provides a procedure to forecast the mechanical behavior of material
in extreme conditions where experimental data is simply not available or very difficult to
collect.

* This chapter is the collaborated work by different research groups. The contribution of the author
of this thesis is to determine with accuracy the necessary input material parameters from atomistic
simulations.

Chapter 6

170

A crucial step in this approach is the appropriate selection and modeling of the
unit processes. These models supply the link between the atomic and mesoscale by
identifying and correlating the relevant material properties, susceptible to atomistic
determination such as energy formation for defects, with the corresponding driving
forces. In this case, we specifically consider the following unit processes: double-kink
formation and thermally activated motion of kinks; the close-range interactions between
primary and forest dislocation, leading to the formation of jogs; the percolation motion of
dislocations through a random array of forest dislocations introducing short-range
obstacles of different strengths; dislocation multiplication due to breeding by double
cross-slip and dislocation pair-annihilation.

A set of material parameters is then obtained from the modeling and identification
stage, which is required to quantify the contribution of each of the unit processes. We
compute these materials properties using a combination of ab initio quantum mechanics
(QM) and force field (FF) calculations. QM describes the atomic interactions from first

principles, i.e., with no input from experiments; unfortunately, QM methods are
computationally intensive and restricted to small systems, making QM calculations
impractical to study most of the materials properties governing plasticity. Force fields
calculations give the total energy of a system as a potential energy function of the atomic
positions and with molecular dynamics (MD) allows the simulation of systems containing
millions of atoms. We used ab initio quantum mechanical calculations (equations of state
of various crystalline phases, elastic constants, energetics of defects, etc.) to develop a
many body force field (FF) (named qEAM FF) for Tantalum. Then, we use the qEAM FF
with MD to calculate the core energy of the l/2a screw dislocation, that of the edge

Chapter6

171

dislocation with Burgers vector b=l/2a in (110) planes. We have also calculated
the formation energies and nucleation lengths of the kinks in b=l/2a screw
dislocations.

One of the appealing features of the present approach is the ability to incorporate
additional unit mechanisms as they may be required by the physics of the problem. For
example, the formation and evolution of dislocation structures are of particular interest in
ductile crystals subjected to large and cyclic deformation. In recent studies, unitmechanism-based micromechanical models have been proposed to elucidate the effective
behavior of dislocation structures on the macroscopic response.

6.2 Unit processes

Plastic deformation in metallic systems is the macroscopic manifestation of
dislocation activity. The resistance to the dislocation motion, therefore, engenders the
hardening properties observed in this type of materials. It is then the complex interplay of
microscopic mechanisms controlling dislocation mobility, dislocation interaction and

dislocation evolution, which confers the macroscopic constitutive properties. In the
present approach, these controlling processes are considered to be orthogonal in the sense
that they are weakly coupled with each other. The interaction among them is only
established through the uniqueness of the macroscopic driving force that is shared, via the
localization process, by all the unit processes.

In this section, we introduce the set of controlling unit processes, which have been
identified for describing the mechanical response of high-purity BCC single crystals, in

Chapter6

172

particular, for Tantalum. We also provide the final expression resulting from the
modeling of each of these processes. A detailed description of the model, including
comparison with experimental data, is given in Ref. 7.

6.2.1 Dislocation mobility: double-kink formation and thermally activated motion of
kinks
We consider the thermally activated motion of dislocations within an obstaclefree slip plane. Under these conditions, the motion of dislocations is driven by an applied

resolved shear stress 't" and is hindered by the lattice resistance, which is weak enough
that it may be overcome by thermal activation. The lattice resistance is presumed to be
well described by a Peierls energy function, which assigns an energy per unit length to
dislocation segments as a function of their position on the slip plane.

In bee crystals, the core of screw dislocation segments relaxes into low-energy
non-planar configurations 5 •8- 14 . This introduces deep valleys into the Peierls energy
function aligned with the Burgers vector directions and possessing the periodicity of the
lattice. At low temperatures, the dislocations tend to adopt low-energy configurations
and, consequently, the dislocation population predominantly consists of long screw
segments. In order to move a screw segment normal to itself, the dislocation core must
first be constricted, which requires a substantial supply of energy. Thus, the energy
barrier for the motion of screw segments, and the corresponding Peierls stress, may be
expected to be large, and the energy barrier for the motion of edge segments to be
comparatively smaller. For instance, Duesbery and Xu 15 have calculated the Peierls stress
for a rigid screw dislocation in Mo to be 0.022 µ, where µ is the <111> shear modulus,

Chapter6

173

whereas the corresponding Peierls stress for a rigid edge dislocation is 0.006 µ, or about
one fourth of the screw value. This suggests that the rate-limiting mechanism for
dislocation motion is the thermally activated motion of kinks along screw segments 16 - 18 .

(a)

(b)

Figure 6-1. Schematic of the double-kink mechanism.

Chapter6

174

At sufficiently high temperatures and under the application of a resolved shear
stress '! > 0, a double-kink may be nucleated with the assistance of thermal
activation 5 •19•20 , and the subsequent motion of the kinks causes the screw segment to
effectively move forward, as shown in Figure 6-1. Under this condition the following
expression for the effective temperature and strain-rate dependent Peierls '!pis obtained:

'f O

f]E kink a sinh( r kink e

/JEkink

(1)

where the effective Peierls stress is given by

£kink

=--bLkink [

(2)

and the reference strain is defined as

r.0kink = 2bn[ P V D'
f-''

(3)

In the preceding equations, b is the Burgers vector, p is the dislocation density, B =
1/kBT, kB is the Boltzmann constant, T is the absolute temperature, and v 0 is the attempt
frequency, which may be identified with the Debye frequency to a first approximation.
Also, [p is the distance between two consecutive Peierls valleys. For bee crystals, [p =

2/3 a if the slip plane is { 110}, [p = ✓
2 a, if the slip plane is { 112}, and [p = ✓
813 a if
the slip plane is {123}, where a is the cubic lattice size 21 . Finally, Ekink is the energy of
formation of a kink-pair and Lkink is the length of an incipient double kink. The formation
energy Ekink and the length Lkin\ which cannot be reliably estimated from elasticity since

175

Chapter6

the energy is composed mostly of core region, can, however, be accurately computed by
recourse to atomistic models as shown in Chapter 5. Modeling of this first unit process
renders the first 2 material properties amenable of atomistic calculations.

0.9
0.8
0.7
0.6

...___ 0.5

0..

0.4
0.3
0.2
0.1

0.05

0.1

0.15

0.2

T kiEkink
Figure 6-2. Temperature dependence of the effective Peierls stress for various strain rates.
Note that the typical order of magnitude of

rtk =10- s6

In Figure 6-2 the dependence of the effective Peierls stress on temperature and
rate of deformation is illustrated. The Peierls stress decreases ostensibly linearly up to a
critical temperature Tc, beyond which it tends to zero. These trends are in agreement with
the experimental observations of Wasserbach 22 and Lachenmann and Schultz23 The
critical temperature Tc increases with the strain rate. In particular, in this model the effect
of increasing (decreasing) the strain rate has an analogous effect to decreasing

Chapter6

176

(increasing) the temperature, and vice versa, as noted by Tang et al.

. In the regime of

very high strain rates ( y > 105 s- 1), effects, such as electron and phonon drag, become
. of ct·IS 1ocat10ns
. 25 ' 26 .
Important
an d contra1 th e ve 1ocity

6.2.2 Dislocation interactions: obstacle-pair strength and obstacle strength

In the forest-dislocation theory of hardening, movmg dislocations could be
impeded by the secondary or "forest" dislocations in their slip planes. As the moving and
forest dislocations intersect, they form jogs or junctions of varying strengths 4 •27 -34 which,
provided the junction is sufficiently short, may be idealized as point obstacles. Moving
dislocations are pinned down by the forest dislocations and require a certain elevation of
the applied resolved shear stress in order to bow out and bypass the pinning obstacles.
For the case of infinitely strong obstacles, the resistance of the forest is provided by the
strength of the obstacle pairs. This obstacle pair strength is subsequently deduced by
considering that point obstacles composing the pair can only provide a finite strength.
The processes imparting the pair-obstacle strength and obstacle strength are described
next.

A. Obstacle-pair strength

We begin by treating the case of infinitely strong obstacles. In this case, pairs of
obstacles pin down dislocation segments, which require a certain threshold resolved shear
stress s in order to overcome the obstacle pair. The lowest-energy configuration of
unstressed dislocation segments spanning an obstacle pair is a step of the form shown as
the thin line in Figure 6-3. Under these conditions, the bow-out mechanism by which a

Chapter6

177

dislocation segment bypasses an obstacle pair may be expected to result m the
configuration shown in Figure 6-3 (bold line).

das
obstacle

Figure 6-3. Bow-out mechanism for a dislocation segment bypassing an obstacle pair.

If the edge-segment length is le, a displacement dae of the dislocation requires a

supply of energy equal to 2 uscrew dae + b r/dge le dae in order to overcome the Peierls
resistance r/dge and to extend the screw segments. The corresponding energy release is
br le dae, Similar contributions result from a displacement das of the screw-segment of

length ls. Retaining dominant terms the obstacle-pair strength is

2Uedge
= 7 screw + __

bl s

(4)

Chapter6

178

The obstacle-pair strength can be therefore estimated by quantifying rp, ls and uectge_ An
expression for the Peierls stress 'rp is given in Eq. (1). The distance between obstacles
along the screw direction l.1 is estimated by statistics assuming a random obstacle
distribution and the core energy per unit length in the edge direction uectge is obtained by
atomistic calculations presented in the following sections.

B. Obstacle strength

In this section we proceed to estimate the obstacle strength that reduces the
obstacle-pair strength described in the previous section. The interaction between primary
and secondary dislocations may result in a variety of reaction products, including jogs
and junctions 4 ·24 •27 -34 _ Experimental estimates of junction strengths have been given by
Franciosi and Zaoui 35 for the twelve slip systems belonging to the family of { 111} planes
and [110] directions in fee crystals, and by Franciosi3 6 for the twenty-four systems of
types {211} [111] and {110}[111] in bee crystals.

The strength of some of these

interactions has recently been computed using atomistic and continuum models 4 •27-29 _
Tang et al. have numerically estimated the average strength of dislocation junctions for
Nb and Ta crystals 24 .

For purposes of the present theory, we specifically concern ourselves with shortrange interactions between dislocations that can be idealized as point defects. For
simplicity, we consider the case in which each intersecting dislocation acquires a jog. The
energy of a pair of crossing dislocations is schematically shown in Figure 6-4 as a
function of some convenient reaction coordinate, such as the distance between the
dislocations. The interaction may be repulsive, resulting in an energy barrier, or

179

Chapter6

attractive, resulting in a binding energy (see Figure 6-4). In the spirit of an equilibrium
theory, here we consider only the final reaction product, corresponding to a pair of jogged
dislocations at infinite distance from each other, and neglect the intermediate states along
the reaction path. In addition, we deduce the strength of the obstacles directly from the
energy supply required to attain the final state, i.e., the jog-formation energy. Despite the
sweeping nature of these assumptions, the predicted saturation strengths in multiple slips
are in good agreement with experiment (cf. Section 6.4), which lends some empirical
support to the theory.

After intersection
Unfavorable Junction _,,/ \.,

_/,i

"--~

" \

Before intersection

.,i

\\ //

Energy inc::re.::ise due
to jog form.::ition

\../ F.::ivorabl e Junction

Figure 6-4. Schematic of energy variation as a function of a reaction coordinate during
dislocation intersection and crossing.

We estimate the jog formation energy as follows. Based on energy and mobility
considerations already discussed, we may expect the preponderance of forest dislocations
to be of screw character, and the mobile dislocation segments to be predominantly of
edge character. We therefore restrict our analysis to intersections between screw and edge
segments. The geometry of the crossing process is schematically shown in Figure 6-5.

180

Chapter6

screw segment (b )
a.

edge segment (b )

screw segment (b )

Before Intersection

After Intersection

Figure 6-5. Schematic of jog formation during dislocation intersection.

Each dislocation acqmres a jog equal to the Burgers vector of the remaining
dislocation. The energy expended in the formation of the jogs may be estimated as
E/xt = bUscrew[l- rcoseaJJ]
(5)

where r = uectge/uscrew is the ratio of screw to edge dislocation line energies. This ratio is
computed by atomistic calculations presented in the next section, renders a value of r
=1.77 for Ta. The resulting jog formation energies for the complete collection of pairs of
{211} and {110} dislocations are tabulated in Table 6-1.

181

Chapter6

Table 6-1. Normalized jog-formation energies resulting from crossings of bee
dislocations.

A2
A.2'
A3
A3'
A6
A6'
B2
B2''
B4
B4'
BS
BS'
Cl
Cl'
C3
C3"

cs
cs··
D1
D1''
D4
04••
D6
1)6''

A2 A2'A3A3'A6A6'B2B2''B4 B4'BSBS'Cl Cl'C3 C3''CS CS'Dl Dl'D4D4'D6 06'
-1.0 1.01.0 1.01.0 1515 15 LS 1.51.5 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4
lO - 1.01.0 1.01.0 3.2 3.2 1212 3.2 3.2 1.81.8 1.81.8 1.81.8 1.81.8 1.81.8 1.81.8
1.0 lO --1.0 1.0 1.0 2.4 2.4 2.4 2.4 2.4 2.4 151.5 1.51.5 1.51.5 2.4 2.4 2.4 2.4 2.4 2.4
1.0 lO 1.0-- 1.0 1.0 1.81.8 1.81. 81.81.8 12 3.2 3.2 3.2 3.2 3.2 1.81.8 1.81.8 1.81.8
1.0lO1.01.0--1.0 2.42.4 2.4 2.4 2.42.4 2.4 2.4 2.4 2.4 2.4 2.4 1.51.5 1.51.5 1.51.5
lO 1.01.01.0 1.0--1.81.8 1.81.81.81.81.81.8 1.81.8 1.81.8 3.2 3.2 3.2 3.2 3.2 3.2
1.5 1.5 151.5 1515 --1.0 lO 1.01.01.0 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4
12 12 3.23.2 3.23.2 1.0-- 1.01.01.01.0 1.81.8 1.81.8 1.81.8 1.81.8 1.81.8 1.81.8
2.4 2.4 2.4 2.4 2.4 2.4 1.01.0 -1.01.01.0 2.4 2.4 2.4 2.4 2.4 2.4 1515 1515 1.515
1.8 1.8 1.81.8 1.81.8 1.01.0 1.0--1.01.0 1.81.8 1.81.8 1.81.8 3.23.2 3.23.2 3.23.2
2.4 2.4 2.4 2.4 2.4 2.4 1.01.0 lO 1.0--1.0 151.5 1.515 1.51.5 2.4 2.4 2.4 2.4 2.4 2.4
1.81.81.81.8 1.81.8 1.01.0 lO 1.01.0--123.2 3.23.2 3.2 3.2 1.81.8 1.81.8 1.81.8
1.8 1.8 1.81.8 1.81.8 1.81.8 1.81. 81.81.8 -1.0 1.01.0 1.01.0 3.23.2 3.23.2 3.23.2
1.8 1.8 1.81.8 1.81.8 1.81.8 1.81. 81.81.8 10-- 1.0 1.0 1.01.0 3.2 3.2 3.2 3.2 3.2 3.2
1.5 1.5 1.51.5 1.51.5 2.4 2.4 2.4 2.4 2.4 2.4 lO 1.0 --1.0 1.01.0 2.4 2.4 2.4 2.4 2.4 2.4
12 12 3.2 3.2 3.2 3.21.81.81.81.81.81.8101.01.0--1.01.0 1.81.8 1.81.8 1.81.8
2.4 2.4 2.4 2.4 2.4 2.4 1.51.5 151.51.51.5 lO 1.0 1.01.0 --1.0 2.4 2.4 2.4 2.4 2.4 2.4
1.81.81.81.8 1.81.8 3.23.2 1212 3.23.2 lO 1.01.01.0 1.0-- 1.81.8 1.81.8 1.81.8
1.81.8 1.81.8 1.81.81.81.8 1.81.81.81.8 3.2 3.2 3.2 3.2 3.2 3.2 --1.0 1.01.0 1.01.0
1.8 1.8 1.8 1.8 1.81.8 1.81.8 1.8 1. 8 1.81.8 3.2 3.2 3.2 3.2 3.2 3.2 1.0-- 1.0 1.0 1.0 1.0
2.4 2.4 2.4 2.4 2.4 2.4 1515 151.51.51.5 2.4 2.4 2.4 2.4 2.4 2.4 1.01.0 --1.0 1.0 1.0
1.8 1.8 1.81.8 1.81.8 3.23.2 1212 3.2 3.2 1.81.8 1.81.8 1.81.8 1.01.0 1.0-- 1.01.0
1.5 1.5 1.51.5 1.5 1.5 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 2.4 1.0 1.0 1.0 1.0 --1.0
12 12 3.23.2 3.2 3.2 1.81.8 1.81. 81.81.8 1.81.8 1.81.8 1.81.8 1.01.0 1.01.0 1.0--

A derivation entirely analogous to that leading to Eq. (1) yields the following
expression for the strength of an obstacle in the slip system a produced by a forest
segment in the system B.

sa/J

r. a

. h( --e (/JE/x'Jl))
-R£iog
-o - a Slll
·a
JJs afJ
Yo

(6)

where the strength at zero temperature is given by
serf]

EJng
af]

bl a Ljunrt '

and the reference strain rate by

(7)

Chapter6

182

(8)

The lengths

za and Ifm,t describe the geometry of the junction as illustrated in Figure 6-

6. These values, which have been estimated to be of the order of few b in the present
case, can also be obtained by atomistic models.

~unct

Figure 6-6. Schematic of a dislocation line overcoming a junction.

6.2.3 Dislocation evolution: multiplication and attrition

The density of forest obstacles depends directly on the dislocation densities in all
slip systems of the crystal. Therefore, in order to close the model, we require an equation
of evolution for the dislocation densities. Processes resulting in changes in dislocation
density include production by fixed sources, such as Frank-Read sources, breeding by
double cross slip and pair annihilation (see Ref. 37 for review; see also Ref. 38, 39, 40,
41, 42, 43). Although the operation of fixed Frank-Read sources is quickly eclipsed by
production due to cross slip at finite temperatures, it is an important mechanism at low

183

Chapter6

temperatures. The double cross slip, fixed Frank-Read sources and pa1r annihilation
mechanisms are considered next.

A. Dislocation multiplication: fixed Frank-Reed and breeding by cross glide

The rate of dislocation multiplication in a given slip system a produced by fixed
Frank-Reed sources and by breeding by cross glide is written as
(9)

where Ao is a constant associated with the fixed Frank-Read production; this parameter is
mere topological than material dependent.

B. Attrition: pair annihilation

The rate of dislocation attrition due to pair annihilation may finally be estimated
as
(10)

where K is the effective annihilation distance. This is the maximum distance at which two
screw segments with opposite direction and forced to move with a velocity v = y I hp
will annihilate. This distance can be estimated by simply equating the time required for
trapping and escaping. Trapping is governed by the elastic interaction forces (attraction)
while escaping by the applied strain rate. Then,

Ka(A+ ✓A2+l),

- = - +---,====where

(11)

184

Chapter 6

A -_ e -{JE

/3'£ jog Yo· jog / Y· a ,

(12)

is a factor depending on the strain rate and temperature,
Y.0;og = 2bn[
f-" P D'

(13)

is a reference slip-strain rate and Kc is the cutoff value corresponding to the effective
screening distance. It follows that the critical pair-annihilation distance K decreases with
increasing strain rate and decreasing temperature. Thus, at high strain rates the
dislocation velocities are high and the probability of being captured by another
dislocation diminishes accordingly. Additionally, an increase in temperature increases the
dislocation mobility and speeds up the annihilation process, which results in an attendant
increase in annihilation rates. The rate of annihilation is then modulated by the nucleation
energy of a jog Ejog, which can be calculated from atomistic simulations.

6.3. Atomistic modeling of dislocations properties

In the previous section, we have identified the following set of material

parameters required to estimate the contribution of each of the controlling unit processes:
. fl y d escn'be th e computation
E kink , Lkink , uedge/Uscrew , an d Ejog . I n t h'1s sect10n,
we b ne
of

these parameters using a first principles based force field with molecular dynamics.

Quantum mechanics (QM) describes the atomic interactions from first principles,
1.e., using no empirical input. Unfortunately QM methods are computationally too
intensive and thus only applicable to small systems (hundreds of atoms) and short times
(picoseconds). The studies of most of the unit processes that govern the plasticity of
materials (such as dislocation mobility, kink energies, etc.) involve many atoms and long

Chapter6

185

simulation times. Such problems require the use of force fields, where the internal energy
of the system is given by a potential energy function of the atomic positions and does not
involve the solution of Schrodinger's equation. The drawback of using potentials to
describe the atomic interactions is that some accuracy is lost; it is thus of great
importance to use accurate force fields to describe the atomic interactions.

We developed a many body force field for Tantalum based on accurate QM
calculations that can be used with molecular dynamics (MD) to simulate systems
containing millions of atoms. We fitted an embedded atom model type force field (named
qEAM FF) to a variety of ab initio calculations, including the zero temperature equation
of state (EOS) for bee, fee, and A15 phases of Ta in a wide pressure range, elastic
constants, vacancy formation energy and energetics of a shear transformation in the
twinning direction. Ta is a bee metal and no phase transition to other crystalline phase is
known, but using QM we can calculate the EOS of thermodynamically unstable or
metastable phases (such as A15, fee, hep, etc.). Including data about these high-energy
phases, with different coordination numbers, in the force field training set is important to
correctly describe the atomic interactions near defects, such as dislocations, grain
boundaries, etc.
We have used the qEAM with MD to study a variety of materials properties44 . We
have calculated the melting curve of Ta in a wide pressure range; the calculated zero
pressure melting temperature T melt = 3150K is in very good agreement with the
experimental result of 3290K; this is an important validation given the fact that the
qEAM FF is based only on zero temperature ab initio data. The calculated thermal

186

Chapter6

expansion is also in good agreement with experimental results. We have also used the
qEAM FF with MD to study spall failure in Ta at high strain rates

We use the qEAM FF to calculate a variety of dislocation properties 14 , such as
core energies, Peierls stress, kink formation energies. As pointed out in previous sections,
these are the fundamental quantities that govern plasticity in metals. The accuracy of the
materials parameters obtained from these calculations is best assessed by their use in
macroscopic models that can be directly compared with experimental results. These
quantities could not be directly measured experimentally. The best validation of the
accuracy of the atomistic calculations is through their use in macroscopic models that can
be directly compared with experimental results. The following subsections describe some
of these calculations; in subsection 6.3.1, we show the calculation of the core energy of
edge and screw dislocations in Ta and in subsection 6.3.2 we calculate the double kink
formation energy and nucleation length.

6.3.1. Core energy of 1/2a screw and edge dislocations

In order to study static properties of the 1/2a screw dislocation in Ta,
such as core structure and energy, we use a dislocation quadrupole in a simulation cell
with periodic boundary conditions. Two of the dislocations have Burgers vector
b=l/2a and the other two have b=-l/2a<-1-1-1>. Such an arrangement of
dislocations minimizes the misfit of atoms on the periodic boundary due to the effects of
periodic images. We build the dislocations using the atomic displacements obtained from
elasticity theory and then we relax the atomic coordinates using the qEAM FF. In the bee
structure, there are two kinds of dislocation core configurations (easy core and hard core)

187

Chapter6

that can be transformed to each other by reversing the Burgers vector. In this work we
focus on the lower energy easy cores. In Figure 6-7 we show the differential
displacement map (DD) of our relaxed quadrupolar system. In the DD maps, atoms are
represented by circles and projected on a (111) plane. The arrows represent the relative
displacement in [111] direction of neighboring atoms due to the dislocation. We can see
from Figure 6-7 that the equilibrium dislocation core obtained using qEAM FF has threefold symmetry and spreads out in three <112> directions on {110} planes.

[111]
[11-2]

O••o•¾O••o••O••O••o••O•~o••
It

•:O •: •

[1-1 OJ

~ ff

O••O••O••O••O••O••O••O••O••

•O••O•-»O••O••O••O••O••O••O

• ~•

1it

* o

• -~o < • -c: • --= o

o • • o • • o • •
• ~. • > o~ • • o • * o

oB
• A

,. ~o ..:: •

O••a••O••O••a••O••O••O••<>••

Figure 6-7. Differential displacement map of a relaxed quadrupole of screw dislocations
in Ta.

188

Chapter 6

Let us define strain energy as the total energy of our system once the perfect
crystal energy is subtracted. The total strain energy of a system containing dislocations
can be divided into two terms: core energy (Ee) and elastic energy (Ee)- The latter
contains the self-energy of each dislocation and their interactions and can be calculated
using linear elasticity theory. The core energy is the energy contained close to the
dislocation line (closer than some distance re called core radius), where, due to the large
strains, elasticity theory is not valid and the details of the interatomic interactions are
important. For our quadrupole system the total strain energy takes the form 13
(14)

where K depends on the elastic constants, d 1 and d2 are the nearest separation of
dislocations along <11-2> and <1-10> directions and A(d 1/d 2) is a geometric factor which
comes from the dislocation interactions.

We studied quadrupolar dislocation cells of different sizes. In Figure 6-8 we
show the minimized energy as a function of ln(d 1/rc)+A(d 1/d 2 ) for the different simulation
cells. We took the core radius to be rc=2.287b; this is a typical value used in previous
studies 11 •13 . We can see from Figure 6-8 that the total energies follow a straight line as
predicted by elasticity theory [Eq. (14)], showing that the value chosen for the core radius
is large enough to take account for the non elastic region near the dislocation line.

From a linear fit to our data we determine the core energy Ec=l.404 eV/b and K =
3.3497 x 10- 2 eV/A. 3 . The value of K can also be computed from the elastic constants
giving 3.3492 x 10- 2 eV/A. 3 in excellent agreement with the one obtained from the fit.
Recent ab initio calculations of core energy (using periodic cells containing 90 atoms)

189

Chapter6

give 0.86 eV/b, lower than the value obtained with qEAM FF and the dislocation cores
· 13 .
are compact an d symmetnc

~2.2

dislocation quadrupole in Ta
qEAM FF

(lJ

'-../

2.1

I-<
(lJ

p..

b()
I-<
(lJ

c::
(lJ

·§ 1.9
l:l
rll

N=2310

B 1.8

1.7
0.6

N=l890

0.7

0.8

0.9

1.1

1.2

1.3

ln(d/2b) + A

Figure 6-8. Total strain energy of the quadrupolar system as a function of
In(~)+ A(~); the number of atoms in each simulation is shown. The line is the linear
r,_.
dz
fit to our atomistic data.

Using the qEAM we can calculate the strain energy associated with each atom. In
Figure 6-9 we show the atomic energy distribution (number of atoms per dislocation per
Burgers vector as a function of their strain energy) for a system containing 5670 atoms in
the periodic cell. We can see that there are 6 atoms with atomic strain energy higher than
0.15 eV and another 6 atoms with energy in the range 0.06-0.08 eV. They correspond to

190

Chapter6

the 12 atoms per dislocation per Burgers vector closer to the dislocation line and their
total energy is 1.400 eV/b, very similar to the core energy obtained from Eq. (14). The
rest of the atoms have lower strain energy and can be considered as the elastic part of the
system. We can then define the dislocation core as formed by the 12 atoms per Burgers
vector with higher energy.

!Quadrupole, N=5670 I
elastic region

core region ( 12 atoms I

170.5

6 atoms

0.()1 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 O.l.HJ.14 0.15 0.16 0.17

straln energy (cV)

Figure 6-9. Histogram of atomistic strain energy distribution for the quadrupolar
arrangement of screw dislocations. The cell contains 5670 atoms and is 7 Burgers vectors
long.

We have also calculated the core energy of the edge dislocation with
b=l/2a on a (110) plane. We build a simulation cell with axis oriented along
<112> (x axis), <110> (y axis), and 1/2a (z axis); this cell contains 6 atoms. We

191

Chapter6

then replicate the cell 3 times along X, 16 times along Y, and 20 times along Z; the
number of atoms in the cell is then N=5760. We then remove 108 atoms to form a dipole
of edge dislocations. Once the system 1s relaxed (both atoms and cell parameters), we
have a 24.3967 A x 75.1824 A x 56.632 A cell. Figure 6-10 shows a snapshot of the
atoms projected on a <112> plane.

Edge dislocations in Ta

.... L · • • LL L . . . . ,. L 0.. ♦♦ .~:.I, ....... 1.•. •••.,_ "• ,_ ......... ~ ...._ .... 11, :_ ._ . . . . .._ 1..1. . . . . .._._,I.. •••

..._, . _ •••.• .._.._.,,..,,._, ~ ....... ..,i...••••.tt •••±• ■ t,t 1o.lll ■■ i

. . . . ! . ._ ~

..... L !

< ,a ■ .. . ._ '- l

- ,. 1-. ■•-.. t. L ■ .... ~ L-.. ■- ■ L -

. . . . . •- I,. l

t...

~, . . . ,i.~,\.•I+-. .,,,

....... ! ... •- ...... L t <.. ■ ... , " ! , I .... ■ L ,.

••+

....

■ W.., L. • - . • . . . . . . •~.._Le ■ • .._ L ... ■• ■ L ~.I, ■■••., t. L .. ■ •

.... L . . . . ._ ""l • ■•• t.L ....... I,.._, ..... ;..,i, . . . . . . . . . . . :..'-,ti.,,. l,.l • ■ •:.1 • .., ■■ •t.t,.., ■••
t...•·•L....._..._ . . . . . . ..._L . . . . L,_L., ■ •!_o,. . . . . . -11 .......... • ...... ~,........ flL>...
'.I, . . . . L ~•,,.•lit~ I,,,,
,.;.I,

•• ■ t..l-·... • ■ •t.

••tt~.

L •-.Ill_, .. L l

•.. :. ■••,._ ..,_ ft•I.
l,.

,. _l·-·••-LL•••c l ....

11,

L···L•.

._aill,..

'11 +9 +If':.,. 'k-L- • ■ • 1.. '- .._ f •• .._ :_ 1,. -t-llJ. L .._ \.. • ·

tllllli+ .

,,:.1.:..• ■ •1. ·•·•• ■ :I,., .... ■■ •..,.-...~a-••

' ~ . . . . '.L.L.•4'·'k .. .._ ........ ~.- . ._.,l.,l,.!

..,_ ........ ,.,_._ . . . . :..LL•4•-..:..~a-•6,,; ■■ -t-0.'>. ■■ ,1

I ..... ...__~ ... •-.._ , ......

••• ■ !-IL. • • • -.. .. :.••••·•·-•··

Figure 6-10. Snapshot of the relaxed edge dipole configuration. The cell contains 5652
atoms.

192

Chapter6

100 -.----......---.---,-----.---,-----.--~----,,--r---,---.----,,------r--,---,

TaqEAMFF
edge dislocation

tl.l

......

E 40
!:l

0 '---'--------------.......,-------L.A...11..-..------.___----. . . . . .______,....-..,...,___.IJ...____.____--'-___,_______.____,

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

strain energy (e V)
Figure 6-11. Histogram of atomistic strain energy distribution for the dipole of edge
dislocations. The number of atoms is given by per dislocation per 1/2a length.

In Figure 6-11 we show the energy distribution for the edge dislocation (number
of atoms per dislocation and per a length as a function of their energy). Figure 6-11
shows that the core of the edge dislocation contains atoms with higher energies and a
broader distribution of energies as compared with the screw case (Figure 6-9). Taking
into account Figure 6-11, we define the core of the edge dislocation as formed by those
atoms with strain energy higher than 0.1 eV. This definition leads to 36 atoms per
a<112> or - 4.42 atoms per A and to a core energy of

E;:;,e= 0.860 eV/A (in the case of

Chapter6

193

the screw we had 12 atoms/b or~ 4.17 atoms per A). The ratio between the core energy
of the edge and that of the screw is E;:;; I E;~;;w = 1.77. It is important to mention that
changing the number of atoms considered to belong to the core changes the core energy,
but the difference is minor. Had we taken the 34 atoms per a with higher energy as
the core (leading to ~ 4.18 atoms I A, a density very similar to the one obtained in the
screw dislocation), we would have gotten the core energy E;,,~:ee= 0.84 eV/ A.

6.3.2. Kink pair energy and nucleation length

As already explained, the kink pair mechanism controls the mobility of screw
dislocations in bee metals and atomistic simulations can provide the details of this
mechanism.

As we can see from Figure 6-7, the core of the screw dislocation spreads in three
<112> directions, this leads to two distinct, but energetically equivalent, core
configurations; we name them as positive (P) and negative (N) cores. The shortest (and
lowest energy) kinks possible involve the displacement of the position of the dislocation
line in the (111) plane from one equilibrium position to a nearest neighbor equilibrium
position; the displacement involved is 1/3 a. There are six possible <112>
directions but only two need to be considered by symmetry, this leads to two kink
directions, which we call left (L) and right (R). The two dislocation cores (N and P) and
two directions (L and R) lead to 8 different single kinks: NRP, NRN, PRP, PRN, NLP,
NLN, PLP and PLN. We have studied all of them in detail; here we will concentrate on
the single kinks that lead to the lowest energy kink pair. We calculated the formation

Chapter 6

194

energy and length of the various kinks using quadrupole arrangements of dislocations as
explained in Chapter 5. The simulation cell lengths are 40.7 A in the [11-2] direction,
42.3 A in the [l-10] direction and 431.8 A in the [111] direction. The whole simulation
cell contains 40,500 atoms. The details of these calculations can be found in Ref. 46. We
calculate the kink energy as the difference of strain energy between the quadrupolar
systems containing kinks and perfect straight dislocations. The energy difference divided
by four is the formation energy for each kink. Using the qEAM FF, we find that the
lowest energy kink pair is formed combining the PLN and NRP kinks. We define the
kink pair nucleation energy as the sum of the formation energy of the two single kinks
leading to Ekink = 0.730 eV. This result is comparable to that obtained by Yang et al.
(0.96 eV) using the quantum-based multi-ion interatomic potentials derived from the
model generalized pseudopotential theory (MGPT). The nucleation energy calculated in
this way does not take into account the attractive interaction between the two kinks that
lowers the nucleation energy. This interaction energy is very small (- 2%) for separation
of kinks larger than - 15 b 11 • 20 .

As explained above, a critical parameter for the micromechanical modeling of
plasticity is, apart from the kink pair energy, its nucleation length Lkink• We studied both
the energetics and structure of the various kinks along the dislocation line. Figure 6-12
shows the extent of the kinks both from structural and energetic points of view. We show
the position of the dislocation in the direction of the kink along the dislocation line for a
PLN kink [Figure 6-12(a)] and NRP kink [Figure 6-12(c)]. We also show the total strain
energy of the quadrupolar system along the dislocation line for the PLN [Figure 6-12(b)]
and NRP [Figure 6-12(d)] kinks. It is calculated by summing the atomic strain energies

195

Chapter6

for atoms in every 1 b slice in the [111] direction. The structural length of the PLN kinks
is L:~N = 8 b [Figure 6-12(a)]; while its "energetic extent" is L:;; = 14 b [Figure 612(b)]. For NRP kinks, we obtain L~:P = 8 b [Figure 6-12(c)] and L::P = 20 b [Figure 612(d)].

Going back to the definitions of the parameters entering the equation that governs
the dislocation mobility [Eq. (1) and Eq. (2)]; the effective Peierls stress (to) in Eq. (2) is
defined as the applied stress for which the nucleation free energy for a kink pair (~G) is
zero. ~G is given by
(15)
where Lkink is the effective kink pair nucleation length and [pis the distance advanced by
the dislocations; in the kinks studied here, [p = I 1/3 a<112> 1- The second term in the
right-hand side of Eq. (15) is the work done by the external stress when the kink is
nucleated. Figure 6-13 shows a schematic diagram of a PLN-NRP kink pair. We can see
that the work done by the external stress to nucleate the kink pair can be divided in four
terms:

LPLN

tb[

Lkink

LPLN

= tb[ ( -2!!:.,_ + ene P

LPLN
sir

LNRP
ene

LNRP

sir

NRP

+ L,,·1r )
2'

(16)

where Lkink is the effective kink pair length. In Figure 6-13 we show the four terms in the
right-hand side of Eq. (16). Note that Eq. (16) assumes that the kinks are straight lines
connecting the two equilibrium positions of the dislocation. This way we obtain the
effective kink pair nucleation length Lkink = 17 b.

Chapter 6

196

i".9

Ta qEAl\.IFF
PLNkink

,.8

$ ~5
~ ~o

p:.J

-~ti

i".5

r.n ,.4

'E0 l

E--<

I .3

f-----\
Sb

- - -i,itbont kink

I.I

~ i".6

Pl.Ill klDk

j,....

-1'1thPLNk

~ --

TaqEAl\'IFF
PLN kink energy

14 b

'-~
80

100

llC

100

110

z (b)

z (b)
40

7.9

Ta qEAl\iIFF
35

~Pklnk

7.8

p:.J

-~ti

7.5

r.n 7.4

'E0 l

E--<

1s.;7

7.3

20b

040

7.7

!IIRP klDk

j.... ~o

- - without klDk

~ 7.6

with !11RP klDk

!IIRP klDk enttllY

,,...__

TaqEAMFF

7.~

Z (b)

100

110

100

110

z (b)

Figure 6-12. PLN and NRP kinks in Ta using the qEAM FF. (a) PLN kink: Dislocation
position in the [11-2] direction along the dislocation line; we can see the dislocation
moves from an equilibrium position to the next in a length of 8 Burgers vectors. (b) PLN
kink: total strain energy in the quadrupolar system with four PLN kinks along the
dislocation line. The system is divided in slices with thickness equal to b and the energy
in each region is calculated. (c) NRP kink: Dislocation position in the [11-2] direction
along the dislocation line; we can see the dislocation moves from an equilibrium position

197

Chapter 6

to the next in a length of 8 Burgers vectors. (d) NRP kink: total strain energy in the
quadrupolar system with four PLN kinks along the dislocation line.

3.5

2.5

,-_

..0

2 3
I I

1.5

=0.5
.·;:J0.5

,-;
tl:l

0..

...

o.r.:

PLNkink
14 b

1.5

NRPkink
20b

2.5 0

20
30
40
position along <111> (b)

Figure 6-13. Schematic diagram of a kink pair formed by a NRP and PLN single kinks.
The four terms entering in the work expression [Eq. (16)] are shown in the figure.

The remaining material parameter is the nucleation energy of a jog ei 0 g. In this
work we take Ejog as the PLN-NRP kink pair nucleation energy.

Chapter6

198

6.4 Experiment, validation and prediction
To test the predictive capabilities of the multiscale approach, we first select a set
of material parameters to best fit the experimental results, then we compare these
parameters against the atomistically computed ones, and finally we predict the
macroscopic response using the atomistic parameters. As we shall see, the agreement
between the fitted and computed by atomistics material parameters is remarkable, and the
predicted macroscopic response retains most of the experimental features. These facts
provide confidence in the multiscale modeling approach, indicating that even in the case
that experimental data would not have been available, still the macroscopic behavior
could have been predicted based only on atomistic calculations.

The experiment data correspond to uniaxial tests on Ta single crystals of Mitchell
and Spitzig47 . In these tests, 99.97% pure Ta specimens were loaded in tension along the
[213] crystallographic axis, at various combinations of temperature and strain rate. In
particular we considered temperatures ranging from 296 K to 573 K, and strain rates
ranging from 10- 1 s- 1 to 10-5 s- 1 . The numerical procedure employed for the integration of
the constitutive equations has been described elsewhere48 . The constitutive update is fully
implicit, with the active systems determined iteratively so as to minimize an incremental
work function. All stress-strain curves are reported in terms of nominal stress and
engineering strain.

Two different sets of material properties were used for the numerical simulations.
The first set was obtained by fitting the simulation results to the experimental results.

199

Chapter 6

Table 6-2 identifies the subset of parameters that are also amenable to direct calculation
by atomistic-based methods. The table lists the parameter values obtained by these
methods, as described in Sections 6.3, in parallel with the values obtained by the fitting
approach. Thus, in the second set of properties that were used for numerical simulations,
atomistic-based values replace fit-based values, when available. This is the case for the
edge and screw dislocation self-energies, as well as the kink-pair formation energy and
length. Clearly, those two sets do not differ by much, which strongly support the validity
of the advertised multiscale paradigm. For a complete list of parameters for the model,
the reader should refer to Ref. 7.
Table 6-2. Material parameters for Tantalum.

Parameter

Fitted set

Atomistic set

Ekink (eV)

0.70

0.730

Lkink/b

uedge/µb2 (*)

0.2

0.216

Uedge/Uscrew

1.77 **

1.77

l lb

rJunct/b

Ecross (eV)

0.67

0.730

AFR

2.3

4.5 ***

1250

500 ***

* µ =-C44 +-(Cll -C,2),

** Taken from the atomistic simulations,
*** Not computed by atomistic simulations.

Chapter6

200

120 , - - - - - - - - - - - - - - - - - - - - - - - - - ,
110
100
90

......
ti:!
a.

::::
( I)
(I)

'tii
ti:!

"><

70
60
50
40
30

. 1 o-3 s -1
c:=

20
10
00

0.1

0.2

0.3
Axial strain

0.4

0.5

(a) Experimental data of Mitchell and Spitzig47

120
110
T = 296 K

100
90

......
ti:!
a.

::::
( I)
(I)

'tii
ti:!

"><

80
70
60
T == 373 K

50
40
T = 398 K

c:· = 1 0- s_,

20
T::: 573 K
10

0.1

0.2

0.3
Axial strain

0.4

0.5

(b) Predictions of the model with fitted parameters

201

Chapter6

120
110
100
90

......
ttS

a.
:l:

..... 70
Ill
Ill

!!!
1ii
ttS

">ii
<(

60
T = 373 K

50
40

T =398 K

30
20

T = 573 K

10
00

0.1

0.2

0.3
Axial strain

0.4

0.5

Figure 6-14. Temperature dependence of stress-strain curves for [213] Ta single crystal

202

Chapter 6

120
110
100
90

......
as

:a:
......

C.
(I)
(I)

I!!
"Iii
as

'><

60
50
40
30

0.1

0.2

0.3

0.4

0.5

Axial strain

(a) Experimental data of Mitchell and Spitzig47

120
110
100
90

......
as

:a:
......

C.
(I)
(I)

I!!
"Iii
as
'>ii

60
50
40

T = 373 K

30
20

10
00

0.1

0.2

0.3

0.4

0.5

Axial strain

(b) Predictions of the model with fitted parameters

203

Chapter 6

120
110
100
90

......
l'l:I
C.

:!:
......
UI
UI

!!!
1ii

80
70
60
50

l'l:I

'<(
><

T = 373 K

30
20

ii::::: 10-5

10
00

0.1

0.2

0.3
Axial strain

0.4

0.5

(c) Predictions of the model with atomistic parameters
Figure 6-15. Strain-rate dependence of stress-strain curves for [213] Ta single crystal
(T=373 K).

Figure 6-14 and Figure 6-15 show the predicted and measured stress-strain curves
for a [213] Ta crystal over a range of temperatures and strain rates. One can compare,
from top to bottom: the experimental results, the results obtained after fitting the
parameters, and the results obtained with atomistic-based parameters. It is evident from
these figures that the model, with both sets of parameters, captures salient features of the
behavior of Ta crystals such as: the dependence of the initial yield point on temperature
and strain rate; the presence of a marked stage I of easy glide, specially at low
temperature and high strain rates; the sharp onset of stage II hardening and its tendency to
shift towards lower strains, and eventually disappear, as the temperature increases or the

Chapter 6

204

strain rate increases; the parabolic stage II hardening at low strain rates or high
temperatures; the stage II softening at high strain rates or low temperatures; the trend
towards saturation at high strains; and the temperature and strain-rate dependence of the
saturation stress. Thus, the predictive approach based on atomistic methods clearly shows
its capacity to produce results matching the experimental evidence.
The theory reveals useful insights into the mechanisms underlying the plastic
deformation behaviors of the single Ta crystal. For instance, since during state I the
crystal deforms in single slip and the secondary dislocation densities are low, the Peierls
resistance dominates and the temperature and strain-rate dependency of yield owe mainly
to the thermally activated formation of kinks and crossing of forest dislocations. It is
interesting to note that during this stage the effect of increasing (decreasing) temperature
is similar to the effect of decreasing (increasing) strain rate, as noted by Tang et al. 24 . The
onset of stage II is due to the activation of secondary systems. The rate at which these
secondary systems harden during stage I depends on the rate of dislocation multiplication
in the primary system. This rate is in tum sensitive to the saturation strain y5a1, which
increases with strain rate and decreases with temperature. As a result, the length of the
stage I of hardening is predicted to increase with strain rate and decrease with
temperature, as observed experimentally. Finally, the saturation stress is mainly governed
by the forest hardening mechanism and, in particular, by the strength of the forest
obstacles. This process is less thermally activated than the Peierls stress, since the
corresponding energy barriers are comparatively higher. Consequently, the stress-strain
curves tend to converge in this regime, in keeping with observation.

205

Chapter 6

The apparent softening observed in simulation results at the lowest temperature
(296 K) and the highest strain rate (10- 1 f 1) is actually an effect of the boundary
conditions, allowing some level of rotation of the specimen. Since in those cases, the
material hardening is relatively low (stage I only), this geometrical softening dominates
in the apparent macroscopic behavior. In the other cases, the activation of several systems
at high strains results in a more isotropic deformation, in tum leading to limited rotations.
In order to take the exact experimental boundary conditions into account, a finite element

model of the whole specimen should be used, allowing for a nonhomogeneous
deformation field.

6.5 Conclusion
In this chapter, we present a modeling approach to bridge the atomistic with
macroscopic scales in crystalline materials. The methodology combines identification and
modeling of the controlling unit processes at microscopic level with the direct atomistic
determination of fundamental material properties. These properties are computed using a
many body force field derived from ab initio quantum-mechanical calculations. This
approach is exercised to describe the mechanical response of high-purity Tantalum single
crystals, including the effect of temperature and strain-rate on the hardening rate. The
resulting atomistically informed model is found to capture the following salient features
of the behavior of these crystals.
1. The dependence of the initial yield point on temperature and strain rate,

2. the presence of a marked stage I of easy glide, specially at low temperatures
and high strain rates; The sharp onset of stage II hardening and its tendency to

206

Chapter6

shift towards lower strains, and eventually disappear, as the temperature
increases or the strain rate decreases,
3. the parabolic stage II hardening at low strain rates or high temperatures,
4. the stage II softening at high strain rates or low temperatures,
5. the trend towards saturation at high strains,
6. the temperature and strain-rate dependence of the saturation stress,
7. the orientation dependence of the hardening rate.

6.6 Comment

The reported multiscale approach for modeling plasticity of Tantalum single
crystal was achieved by strong collaborating among different research groups. It should
be mentioned that the mesoscale simulation results were calculated by A.M. Cuitifio
(Department of Mechanical and Aerosapce Engineering, Rutgers University, Piscataway,
NJ 08854, USA), L. Stainier (Laboratoire de Techniques Aeronautiques et Spatiales,
University of Liege, 4000 Liege, Belgium), and M. Ortiz (Graduate Aeronautical
Laboratories, California Institute of Technology, Pasadena, CA 91125, USA). The
contribution of the author of this thesis is to determine with accuracy the necessary input
material parameters from atomistic simulations.

Chapter 6

207

6. 7 References
1. V. V. Bulatov and L. P. Kubin, Dislocation modeling at atomistic and mesoscopic
scales, Current Opinion in Solid State and Materials Science, 3 (6): 558-561, 1998.
2. R. Phillips, Multiscale modeling in the mechanics of materials, Current Opinion in
Solid State and Materials Science, 3 (6): 526-532, 1998.
3. G. H. Campbell, S. M. Foiles, H. C. Huang, D. A. Hughes, W. E. King, D. H. Lassila,
D. J. Nikkel, T. D. de la Rubia, J. Y. Shu, and V. P. Smyshlyaev, Multiscale modeling
of poly-crystal plasticity: a workshop report, Materials Science and Engineering A,
251 (1-2): 1-22, 1998.
4. R. Phillips, D. Rodney, V. Shenoy, E. Tadmor, and M. Ortiz, Hierarchical models of
plasticity: dislocation nucleation and interaction, Modeling and Simulation in
Materials Science and Engineering, 7 (5): 769-780, 1999.
5. J. A. Moriarty, W. Xu, P. Soderlind, J. Belak, L. H. Yang, and J. Zhu, Atomistic
simulations for multi scale modeling in bee metals, Joumal of engineering Materials
and Technology-Transactions of the ASME, 121 (2): 120-125, 1999.
6. M. I. Baskes, The status role of modeling and simulation in materials science and
engineering, Current Opinion in Solid State and Materials Science, 4 (3): 273-277,
1999.
7.

A. M. Cuitino, L. Stainer and M. Ortiz, Micromechanical modeling of hardening,

rate sensitivity and thermal softening in bee single crystals, Journal of the Mechanics
and Physics of Solids, 2001.

Chapter6

208

8. M. S. Duesbery, V. Vitek and Bowen, Proceedings of the royal Society of London,
A332: 85, 1973.
9. V. Vitek, Proceedings of the royal Society of London, A352: 109, 1976.
10. V. Vitek, Structure of dislocation cores in metallic materials and its impact on their
plastic behavior, Progress in Material Science, 36: 1-27, 1992.
11. W. Xu and J. A. Moriarty, Atomistic simulation of ideal shear strength, point defects,
and screw dislocations in bee transition metals: Mo as a prototype, Physical Review B
- Condensed Matter, 54 (10): 6941-6951, 1996.
12. M. S. Duesbery and V. Vitek, Plastic anisotropy in bee transition metals, Acta
Materialia, 46 (5): 1481-1492, 1998.
13. S. Ismail-Beigi and T. A. Arias, Ab initio study of screw dislocations in Mo and Ta: a
new picture of plasticity in bee transition metals, Physical Review Letters, 84 (7):
1499-1502, 2000.
14. G. Wang, A. Strachan, T. Cagin, and W. A. Goddard III, Molecular dynamics
simulations of 1/2a screw dislocation in Ta, Materials Science and Engineering
A, 309-310: 133-137, (2001).

15. M. S. Duesbery and W. Xu, The motion of edge dislocations in body-centered cubic
metals, Scripta Materialia, 39 (3): 283-287, 1998.
16. P. B. Hirsch, In 5 th International Conference on Crystallography, 139, Cambridge
University, 1960.
17. A. Seeger and P. Schiller, Acta Metallurgica, 10: 348, 1962.
18. J.P. Hirth and J. Lathe, Theory of Dislocations, McGraw-Hill, New York, 1968.

Chapter 6

209

19. J. P. Hirth and R. G. Hoagland, Nonlinearities in the static energetics and in the
kinematics of dislocations, PhysicaD, 66 (1-2): 71-77, 1993.
20. W. Xu and J. A. Moriarty, Accurate atomistic simulations of the Peierls barrier and
kink-pair formation energy for <111> screw dislocations in bee Mo, Computational
Materials Science, 9 (3-4): 348-356, 1998.
21. A. Seeger and L. Hollang, The flow-stress asymmetry of ultra-pure molybdenum
single crystals, Materials Transactions JIM, 41 (1): 141-151, 2000.
22. W. Wasserbach, Philosophical Magazine A, 53: 335, 1986.
23. R. Lachenmann and H. Schultz, Scripta Metallurgica, 4: 33, 1970.
24. M. Tang, B. Devincre, and L. P. Kubin, Simulation and modeling of forest hardening
in body center cubic crystals at low temperature, Modeling and Simulation in
Materials Science and Engineering, 7 (5): 893-908, 1999.
25. T. Suzuki, S. Takeuchi, and H. Yoshinaga, Dislocation Dynamics and Plasticity,
Springer-Verlag, 1991.
26. A. D. Brailsford, Electronic component of dislocation drag m metals, Physical
Review, 186: 959-961, 1969.
27. M. I. Baskes, R. G. Hoagland, and T. Tsuji, An atomistic study of the strength of an
extended dislocation barrier, Modeling and Simulation in Materials Science and
Engineering, 6 (1): 9-18, 1998.
28. D. Rodney and R. Phillips, Structure and strength of dislocation junctions: An atomic
level analysis, Physical Review Letters, 82 (8): 1704-1707, 1999.

Chapter 6

210

29. V. B. Shenoy, R. V. Kukta, and R. Phillips, Mesoscopic analysis of structure and
strength of dislocation junctions in fee metals, Physical Review Letters, 84 (7): 14911494, 2000.
30. G. Danna and W. Benoit, Dynamic recovery of the microstructure of screw
dislocations in high-purity bee metals, Materials Science and Engineering A, 164 (12): 191-195, 1993.
31. M. Rhee, H. M. Zbib, J. P. Hirth, H. Huang, and T. D. de la Rubia, Models for
long/short range interactions and cross slip in 3-D dislocation simulation of bee single
crystals, Modeling and Simulation in Materials Science and Engineering, 6 (4): 467492, 1998.
32. H. C. Huang, N. Ghoniem, T. D. de la Rubia, M. Rhee, H. Zbib, and J. Hirth,
Stability of dislocation short-range reactions in bee crystals, Journal of Engineering
Materials and Technology - Transactions of the ASME, 121 (2): 143-150, 1999.

33. L. P. Kubin, B. Devincre, and M. Tang, Journal of Computer Aided Material Design,
5: 31, 1998.
34. H. M. Zbib, T. D. de la Rubia, M. Rhee, and J. P. Hirth, 3-D dislocation dynamics:
stress-strain behavior and hardening mechanisms in fee and bee metals, Journal of
Nuclear Materials, 276: 154-165, 2000.

35. P. Franciosi and A. Zaoui, Miltislip in fee crystals: a theoretical approach compared
with experimental data, Acta Metallurgica, 30: 1627, 1982.
36. P. Franciosi and A. Zaoui, Glide mechanisms in bee crystals: an investigation of the
case of a-iron through multislip and latent hardening tests, Acta Metallurgica, 31:
1331, 1983.

Chapter6

211

37. D. Kuhlmann-Wilsdorf, Theory of plastic deformation: properties of low energy
dislocation structures, Materials Science and Engineering A, 113: 1, 1989.
38. W. G. Johnston and J. J. Gilman, Dislocation velocities, dislocation densities and
plastic flow in lithium fluoride crystals, Journal of Applied Physics, 30: 129, 1959.
39. W. G. Johnston and J. J. Gilman, Dislocation multiplication in lithium fluoride
crystals, Journal of Applied Physics, 31: 632, 1960.
40. P. P. Gillisand J. J. Gilman, Dynamical dislocation theory of crystal plasticity. II.
Easy glide and strain hardening, Journal of Applied Physics, 36: 3380, 1965.
41. U. Essmann and M. Rapp, Slip in copper crystals following weak neutron
bombardment, Acta Metallurgica, 21: 1305, 1973.
42. K. P. D. Lagerlof, On deformation twinning in bee metals, Acta Metallurgica et

Materialia, 41 (7): 2143-2151, 1993.
43. H. Dybiec, Model of the early deformation stage of bee metals in low-temperature,

Zeitschriftfur Metallurgica, 86 (7): 512-517, 1995.
44. A. Strachan, T. Cagin, 0. Gulseren, S. Mukherjee, R. E. Cohen, and W. A. Goddard
III, First principles force field for metallic tantalum, Physical Review B, submitted.
45. A. Strachan, T. Cagin, and W. A. Goddard III, Critical behavior in spallation failure
of metals, Physical Review B, 63: 0601034, 2001.
46. G. Wang, A. Strachan, T. Cagin, and W. A. Goddard III, Kinks in a/2<111> screw
dislocation in Ta, Journal of Computer Aided Materials Design, in press.
47. T. E. Mitchell and W. A. Spitzig, Three-stage hardening in tantalum single crystals,

Acta Metallurgica, 13: 1169-1179, 1965.

Chapter6

212

48. M. Ortiz and L. Stainier, The variational formulation of viscoplastic constitutive
updates, Computer Methods in Applied Mechanics and Engineering, 171 (3-4): 419444, 1999.