Journal of the Earth Simulator, Volume 7, March 2007, 17–21

First-principles Calculation on Screw Dislocation Core Properties in
BCC Molybdenum
Futoshi Shimizu1*, Shigenobu Ogata2, Hajime Kimizuka3,
Takuma Kano1, Ju Li4, and Hideo Kaburaki1
1 Center for Computational Science & E-systems, Japan Atomic Energy Agency, Tokyo 110-0015, Japan
2 Center for Atomic and Molecular Technologies, and Department of Mechanical Engineering,Osaka University, Osaka 565-0871, Japan
3 Nanotechnology Center, Science Department, JRI Solutions, Ltd., Tokyo 102-0074, Japan
4 Department of Materials Science and Engineering, The Ohio State University, Columbus, OH 43210, USA(Received January 16, 2007; Revised manuscript accepted February 28, 2007)
Abstract Predicting atomistic properties of a dislocation is a first step toward an understanding of plastic
behavior of materials, in particular BCC metals. The core structure and Peierls stress of a screw dislocation in
BCC metals have been studied over the years using the first-principles and empirical methods, however, theirconclusions vary due to the inefficiency of the methods. We have executed first-principles calculations basedon the density functional method, employing the most accurate 1 × 1 × 20 k-point samplings, to determine thecore structure and Peierls stress of the a0/2[111] screw dislocation of molybdenum. We have concluded thatthe core has a 6-fold structure, and determined the Peierls stress of 1.8 GPa for the simple shear strain alongthe (–110)<111> direction.
Keywords: first-principles calculation, screw dislocation, Peierls stress, BCC, molybdenum

1. Introduction
Hydrogen embrittlement in BCC metals is one of themost important and fundamental phenomena associatedwith the environmental assisted cracking observed in structural materials. Before studying the interaction of hydrogenatoms and a dislocation, it is first necessary to understand abehavior of a screw dislocation in BCC metals, in particular, the dislocation core structure and Peierls stress.
Numerous studies using first-principles and empiricalatomistic methods have been carried out to determine theaccurate core properties, however, the results of the core
structure and Peierls stress differ depending on the simulation methods. Two core structures have been proposed
for a screw dislocation in BCC molybdenum by different
computational approaches [1][2][3]; the 3-fold structure(Fig. 1(a)) is obtained using the Finnis-Sinclair potential
and the modified generalized pseudopotential theory(MGPT) potential, while the 6-fold structure (Fig. 1(b)) is
obtained by the tight-binding and the density functional
theory (DFT) methods. The estimated value of the Peierls
potential depends on simulation methods, ranging from2 GPa to 4 GPa.

(a) 3-fold structure(b) 6-fold structure
Fig. 1 Differential displacement map of two proposed screw
dislocation core structures in BCC molybdenum.

There are essentially two factors to be paid attention in
determining the core structure of a screw dislocation using
the atomistic method. One is the boundary condition, andthe other is the interatomic potential. Since a screw dislo-* Corresponding author: Futoshi Shimizu, Center for Computational Science & E-systems, Japan Atomic Energy Agency, 6-9-3
Higashi-ueno, Taito-ku, Tokyo 110-0015, Japan. E-mail: shimizu.futoshi@jaea.go.jp

17

First-principles Calculation on Screw Dislocation in Molybdenum

cation has a long-range strain field around the core, thesimulation region should be taken as large as possible toreduce the effects of the boundaries. In order to take alarge simulation region, an empirical potential, such as,
the embedded atom method (EAM), is introduced toreduce the computational burden of calculating the interatomic forces. However, the empirical potentials are basically constructed by fitting the parameters to the equilibrium states, and there is no guarantee for the empiricalpotential to describe the disordered core region accurately.
On the other hand, although the first-principles calculation
based on the density functional theory gives the mostaccurate results even for the core structure, the computational region is limited to a very small area. Woodward
and Rao [2] proposed a boundary condition for the supercell containing one screw dislocation, in which atomic
displacements due to a screw dislocation in the cell areallotted by the Green function method. Recently, Li et al.
[3] introduced a pair of screw dislocations in a skew
supercell in such a way as to cancel the strains among dislocations effectively. In this paper, we perform the DFT
calculations using the latter supercell setup and report theresult of the most accurate determination of Peierls stressfor the screw dislocation in BCC molybdenum.

2. Method
The detailed structure of electronic and atomic bonding
states is required to understand fundamental mechanicalproperties of materials [4][5]. For the electronic simulations, the DFT is the most dependable method to determine the core structure of dislocations. We use theVienna ab initio Simulation Package [6] (VASP), basedon the DFT with the ultrasoft pseudopotential and thegeneralized gradient approximation (GGA), and implement this on the Earth Simulator (ES), which is a highly
parallel vector supercomputer system of the distributedmemory type, and consists of 640 processor nodes connected by 640 × 640 single-stage crossbar switches. The
re-parallelized VASP code in combination with the highperformance both of the processors and the interconnection network of the ES enables the DFT calculation of theunprecedented scale, such as 231 atoms and 2 × 2 × 40 kpoints, as described below.

2.1 Reliability test using the large supercell
There are two conditions to be considered for the accurate calculation of the DFT method. One is the energy
cutoff convergence for the plane wave expansion of the
wavefunction, and the other is the k-point sampling inusing the supercell. For the energy cutoff convergence,
we have adopted the value of 233 eV since this value isconfirmed to be highly accurate for the calculation of

18

Table 1 The total energy for the BCC Mo perfect lattice structure with 250 atoms.
k-points

Total energy [eV]

1×1×1–2695.249

2×2×2–2702.430

3×3×3–2704.241

4×4×4–2704.163

5×5×5–2704.012

6×6×6–2704.055

Energy difference [eV]–7.182–1.8100.078
0.150–0.042elastic constants [5] and the energy cutoff convergence
depends mostly on the atom species. As a preliminary
examination to check the accuracy of calculation for thesystem of Mo, we performed the relaxation of the BCC
Mo perfect lattice structure with 250 atoms by increasing
the k-points. The results are shown in Table 1. Increasing
the k-point samplings to the highly accurate 6 × 6 × 6 kpoints, we have achieved the convergence accuracy of thetotal energy in nearly 0.1 eV in the relaxation process.
With the 6 × 6 × 6 k-points, we have derived the phonon
dispersion curve of BCC Mo, where we have successfully
obtained the detailed structure near the H point which
agrees well with experiment [7] as shown in Fig. 1.

2.2 Re-parallelization of the VASP code
The “hot spot” of the DFT calculation contains a triple
loop for traversing spin states, k-points, and energy bands.
Although the parallelization over k-points is the simplest
and the most effective, the original VASP code offers only
a parallelization program over the energy bands and theplane wave coefficients, where energy bands are distributed among a group of processors. This parallelization
strategy is predicted to be employed because the VASP
code assumes a small number of k-point samplings inusual calculations. To achieve a high efficiency by reducing communication overhead, the residual minimization
method and the direct inversion in the iterative subspace(RMM-DIIS) matrix diagonalizer algorithm [8], which
reduces the number of orthonormalization steps considerably, should be used since it is the only algorithm whichworks for the distributions simultaneously over energy
bands and plane wave coefficients. However, in the case of
the present BCC molybdenum system consisting ofapproximately 200 atoms with a pair of screw dislocations,
it was found that the RMM-DIIS algorithm has a poor convergence and the procedure requires exceedingly vast
amounts of calculation time. Since the present simulation

J. Earth Sim., Vol. 7, Mar. 2007, 17–21

F. Shimizu

Table 2 Parallel performance of the original and k-point parallelized VASP code.
Number of Nodes (NPROC)
Original code (KPAR=1)
Elapsed Time (sec)
Average Vector Length
Vector Operation Ratio (%)
GFLOPS (%peak)

1 (8)2 (16)

5 (40)

8 (64)

10 (80)38828.3
229.3
99.3
42.5(66.4)

20746.0
224.0
99.2
81.2(63.4)

11309.2
212.2
99.0
169.1(52.8)

9533.2
197.7
98.8
240.4(47.0)

10165.1
185.9
98.7
285.9(44.7)

219751.7
229.0
99.3
83.6(65.3)

58138.1229.2
99.2
203.0(63.4)

k-point parallelized: KPAR
Elapsed Time (sec)
Average Vector Length
Vector Operation Ratio (%)
GFLOPS (%peak)Phonon Dispersion of Mo
8

[ξ00]

[ξξξ]

[ξξ0]

7

80

1.0
0.8

70
0.6

60

k-parallel

0.4 Efficiency

50

20

40

60

80

40
30

6
ν [THz]dislocation using 1 × 1 × 20 k-points. It is confirmed thathigh vector operation ratio of 99.3% using a single node
results in an excellent performance of 66.4% of the peak
performance of the ES. We also found that the k-point
parallelization has almost no influence on the average
vector length; the average vector length is 229.3 using asingle node and 229.1 using 10 nodes, where the maximum length of the ES is 256. The total performance using
10 nodes is 392.9 GFLOPS. The advantage of the k-point
parallelized program over the original program is obvious; the former achieves almost a linear speedup while
the original program has reached its plateau using 64
processes as shown in Fig. 3.

Speeduprequires more k-points, we employed the Davidson block
iteration algorithm by re-parallelizing the program in termsof k-point samplings. The range for the outermost loop of
the total number of k-points, NKPTS, is divided and distributed among a group of processes in a block distribution
scheme; each process group is responsible for the part ofcalculation for NKPTS/KPAR k-points, where KPAR is
the number of process groups. The innermost loop, meanwhile, is still performed in parallel by the original scheme.
Furthermore, the initialization code is modified to allocate
only the part of the array variables for wavefunctions toaccomplish the efficient use of distributed memory. In
order to exploit the vector processors on the ES systemeffectively, we have tuned the program by using the onedimensional FFT routine in the numerical library ASL/ES.
Table 2 shows the performance evaluation results of
the original program and k-point parallelized program on
the ES. Here, the DFT calculations are done for thesupercell containing 231 Mo atoms with a pair of screw

104210.1229.1
99.1
392.9(61.4)

original

20

5

10
4

10

20

30

40

50

60

70

80

3

Processes

2

Fig. 3 Performance of k-point parallelized VASP code.
The speedup of (red curve) the original program and (bluecurve)the k-point parallelized program with KPAR = 1, 2,
5, and 10 at NPROC = 8, 16, 40, and 80, respectively. The
calculations are done for the supercell containing 231
molybdenum atoms with a pair of screw dislocation using
1 × 1 × 20 k-points. The data are normalized as 8 at
NPROC (number of processors) = 8. (inset) the parallel
efficiency normalized as 1 at NPROC = 8.

1
0
Γ

H

P

Γ

N

Fig. 2 Phonon dispersion curve of BCC Mo 250 atoms with
6 × 6 × 6 k-points (lines) compared with the experimental data [7] (open circles). The detailed structure near the
H point has been reproduced.

J. Earth Sim., Vol. 7, Mar. 2007, 17– 21

19

First-principles Calculation on Screw Dislocation in Molybdenum

3. Results
We have examined the effects of k-point samplings
using the supercell containing 231 atoms with a pair of the
3-fold structured screw dislocations. The initial setup is
shown in Fig. 4. This system is scrutinized in calculating
the accuracy of the core energy by changing the size of thesupercell, and it is concluded that the size of 231 atoms isgood enough [3]. In Table 3, the results of the relation
between the total energy and the number of k-points areshown. We find that the results with 2 × 2 × 40 k-points
have already converged to the energy difference of 4.7
meV corresponding to 0.02 meV per one atom. Then thedetailed studies on the core structure and the determinationof the Peierls stress using the k-points of 1 × 1 × 20 havebeen performed.

3.1 Core structure
Starting from the initial supercell with a pair of the 3fold core structured screw dislocations, the system isrelaxed under the condition of fixed supercell shape. It isfound that the initial configuration turns into a stable 6fold core structure, and there exists no energy barrier
between the two structures. This indicates that the 6-foldcore structure is stable, and the 3-fold core structure is noteven meta-stable. The energy difference between the 3fold and the 6-fold core structures for the whole system isestimated to be 0.798 eV, which is 0.399 eV for each dislocation. Then, using the k-points of 2 × 2 × 40, we haveevaluated the forces on each atom and the total energy forthe systems containing the initial 3-fold and the 6-foldcore structures, respectively, keeping the positions of theatoms fixed. As a result, in the case of the 3-fold core
structure, the forces on each atom using the 2 × 2 × 40 kpoints correspond to those using 1 × 1 × 20 k-points withthe energy order of 0.01 eV per each atom. The accuracyof the results for the case of the 6-fold core structure
between the two cases of k-points is also within the energy of 0.01 eV per each atom. On the basis of theseresults, we conclude that the 6-fold core structure shoulde2 -k-points

Total energy [eV]

1 × 1 × 20

–2493.898

2 × 2 × 40

–2494.056

3 × 3 × 60

–2494.051

Energy difference [eV]–0.1580.005

be obtained by the relaxation process using the k-pointsof 2 × 2 × 40 and that the core of a screw dislocation in
the BCC molybdenum has a 6-fold structure.

3.2 Peierls stress
Using the initial configuration of the 6-fold core structure obtained in the relaxation simulation, we applied the
additional shear strain along the e3 direction in Fig. 4 tostudy the stress-driven instability of the dislocation core.
The parameter x in Fig. 4 corresponds to the shear strain
γ 23. With each value of x being increased, atoms arerelaxed to a minimum energy state by the conjugate gradient or the quasi Newton's algorithm. Figure 5 shows theresult of calculation for stress-driven instability of thedislocation core in the (–110)<111> slip system of BCC
molybdenum. At a critical value of xc ≈ 0.146, which corresponds to the shear strain of γc ≈ 1.6%, it is found thatthe core structures can no longer be stabilized and thePeierls stress is estimated to be τp ≈ 1.8 GPa. At the valueof x = 0.16, the two cores in the supercell have moved
toward each other by a0/3[11-2], and the shear stress aswell as the energy of the system decreases significantly(Fig. 5 (c)).

4. Discussion
There have been two proposed core structures of ascrew dislocation as a result of the calculations using theempirical potentials, tight binding method, and the density functional theory. The density functional theory is theh2 = 3.5e1 + 5.5e2 + (0.5+x)e3

a 0[110]e3

Table 3 The total energy for the model system with a pair ofscrew dislocations having 3-fold core structure.

e1-

a 0[112]a 0/2[111]h3 = e 3

h1 = 7e1

Fig. 4 The BCC Mo supercell consisting of 231 atoms.

20

J. Earth Sim., Vol. 7, Mar. 2007, 17–21

F. Shimizu

and the lowest τp value of 1.8 GPa in BCC molybdenum.
Although we use the different boundary condition and kpoints sampling, our result generally comes close to the
DFT result of Woodward and Rao [2].

5. Conclusion

Fig. 5 Stress-driven instability of dislocation cores: (a) an initial position, (b) a critical position, and (c) a displaced
position of the cores. The parameter x corresponds to theshear strain γ23. Positions of dislocation cores are seen as
highlighted atoms in the inset images.

most dependable method, but its k-point convergence wasstill unclear. In order to find a definitive solution to thisproblem, the following three points must be considered.
First, we need to have a well-defined boundary condition
for this problem. This has been already checked and thesystem consisting of 231 atom supercell, as shown in Fig.
4, is found to be accurate for this purpose[3]. Secondly,
we need to evaluate the accurate total energy for this system. The VASP code, which implements the density
functional theory with the ultrasoft pseudopotential, suits
this purpose. Thirdly, computational conditions have to
be clarified by the careful error estimation, including theplane wave and k-point convergence checks.
Under the same supercell setup, Li et al. [3] gave thePeierls stress of τp ≈ 2.4 GPa and ≈ 3.8 GPa using theFinis-Sinclair potential and the tight-binding model, respectively. Using the different first-principles Greens function
boundary condition where the supercell contains one dislocation, Xu and Moriarty [1] used the MGPT and obtained
the value of τp ≈ 3.4 GPa, while Woodward and Rao [2]derived τp ≈ 2.1 GPa by the DFT calculation with four special k-points. We do not find any correlation between core
structures and Peierls stresses predicted by various methods, that is, higher and lower τp values are estimated withthe core structure of both the 3-fold and 6-fold.
In this work, we have performed the most accurate
DFT calculations with twenty k-points using the periodic
boundary condition where the supercell contains twoscrew dislocations, and obtained the 6-fold core structure

J. Earth Sim., Vol. 7, Mar. 2007, 17– 21

We have implemented the VASP (Vienna ab-initio
simulation package) code on the Earth Simulator by parallelization in terms of k-point samplings and achieved agood parallel performance. Using this code and thelargest 1 × 1 × 20 k-point samplings, we have performed
DFT calculations and obtained Peierls stress of 1.8 GPafor the a0/2[111] screw dislocation in BCC molybdenum
under the simple shear stress condition in the(–110)<111> slip system.
(This article is reviewed by Dr. Tetsuya Sato.)

References[1] W. Xu, and J. A. Moriarty, Accurate atomistic simulationsof the Peierls barrier and kink-pair formation energy for<111> screw dislocations in bcc Mo, Comput. Mater. Sci.,
vol.9, pp.348–356, 1998.
[2] C. Woodward, and S. I. Rao, Ab-initio simulation of isolated screw dislocations in bcc Mo and Ta, Phil. Mag. A,
vol.81, pp.1305–1316, 2001.
[3] J. Li, C. Wang, J. Chang, W. Cai, V. Bulatov, K. Ho, andS. Yip, Core energy and Peierls stress of a screw dislocation in bcc molybdenum: A periodic-cell tight-binding
study, Phys. Rev. B, vol.70, pp.104113–1–8, 2004.
[4] S. Ogata, J. Li, and S. Yip, Ideal pure shear strength ofaluminum and copper, Science, vol.298, pp.807–811,
2002.
[5] S. Ogata, J. Li, N. Hirosaki, Y. Shibutani, and S. Yip,
Ideal shear strain of metals and ceramics, Phys. Rev. B,
vol.70, pp.104104–1–7, 2004.
[6] G. Kresse, and J. Furthmüller, Efficiency of ab-initio totalenergy calculations for metals and semiconductors using aplane-wave basis set, Comput. Mat. Sci., vol.6, pp.15–50,
1996.
[7] B. M. Powell, P. Martel, A. D. B. Woods, Phonon properties of niobium, molybdenum, and their alloys, Canad. J.
Phys., vol.55, pp.1601–1602, 1977.
[8] G. Kresse, and J. Furthmüller, Efficient iterative schemes
for ab initio total-energy calculations using a plane-wave
basis set, Phys. Rev. B, vol.54, pp.11169–11186, 1996.

21