Development of n-body potentials for hcp–bcc and fcc–bcc binary transition metal systems

Computational Materials Science 43 (2008) 1207–1215

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Development of n-body potentials for hcp–bcc and fcc–bcc binary transition metal systems J.H. Li a, Y. Dai a, X.D. Dai b, T.L. Wang a, B.X. Liu a,* a b

Advanced Materials Laboratory, Department of Materials Science and Engineering, Tsinghua University, Beijing 100084, China Research Institute of Chemical Defense, Beijing 102205, China

a r t i c l e

i n f o

Article history: Received 18 December 2007 Received in revised form 14 March 2008 Accepted 14 March 2008 Available online 25 April 2008 PACS: 34.20.Cf 83.10.Rs 61.43.Dq 02.70.Ns Keywords: n-body potentials Transition metal systems Metallic glass

a b s t r a c t Within the framework of the second-moment approximation of the tight-binding theory, n-body potentials are proposed for hcp, fcc and bcc transition metals and their alloys. Both the energies and their derivatives calculated from the proposed potentials go smoothly to zero at cutoff radii, thus avoiding the unphysical behaviors that may emerge in simulations. With the assistances of ab initio calculations, the proposed potentials are then applied to Zr, Hf, Cu, V, Nb, Ta and their binary alloys. It turns out that the n-body potentials can well predict the energy sequence of stable and metastable structures of these metals. The vacancy formation energies, surface energies and melting points derived from the potentials also match well with experimental results. Based on the constructed potentials, molecular dynamics simulations reveal that Cu–Nb and Zr–Nb metallic glasses could be formed within the composition range of about 15–72 and 8–80 at.% Nb, respectively, matching well with experimental observations. Voronoi analyses reveal that the dominating atomic packings in Cu–Nb metallic glasses are the icositetrahedron (CN = 14), icosihexahedron (CN = 15) and icosidihedron (CN = 13) with fractions of 33%, 26% and 21%, respectively. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction In empirical or semi-empirical calculations and simulations, interatomic potentials play vitally important roles. Whether simulations can reasonably reproduce the experimental observations entirely depends on the relevance of the employed interatomic potentials. At the early stage, two-body potentials, e.g. the Lennard-Jones potential [1], were applied in a wide range of research. Those pairwise potentials do not incorporate the many-body effect, hence have some inherent drawbacks. In the 1980s, significant progress was made with developing n-body potentials, especially for metallic systems, and a wide variety of n-body potentials have been proposed [2]. Important versions among these n-body potentials are the embedded atom method (EAM) and the secondmoment approximation of tight-binding (SMA-TB) [3,4]. These n-body potentials share a principal view that the total energy of an atom is composed of cohesion energy and de-cohesion energy and that the de-cohesion energy can be represented by a pair term to reflect the electrostatic repulsion. The cohesion energy of an atom, e.g. embedded energy in the EAM potential, is a function of the local electron density at the site of the atom, which is approx* Corresponding author. E-mail address: [email protected] (B.X. Liu). 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.03.020

imately considered as the linear superposition of the contributions from its neighboring atoms. In the SMA-TB potential, the cohesion energy is approximately proportional to the square root of the second-moment of the local density of states, which is the sum of the square of the hopping or transfer integrals between the atom and its neighboring atoms [5]. From the empirical viewpoint, there is not essentially difference between the EAM and SMA-TB potentials in nature. The energy of an atom can be considered as being consists of two parts: pair term and many-body term. The many-body term is a function of the electron density. In some cases, the potential functions may be derived rigorously or approximately from physical arguments. However, in many cases they are guessed experimentally and the parameters are determined through fitting to physical properties. According to the forms of their functions, these n-body potentials may be loosely classed into three groups. The first group is known as the EAM potentials, such as the Sandia (National Lab) potentials developed by Daw et al. [6–9], Johnson potential [10,11], Cai and Ye potential [12], Wadley potential [13], Mishin et al. potential [14,15], in which the embedded functions are logarithmic form. The second group is the SMA-TB potentials, such as the Rosato– Guillope–Legrand potential [4], Willaime–Massobrio potential [16] and Cleri–Rosato potential [17], in which the pair and electron density terms are exponential forms. The third group is the

1208

J.H. Li et al. / Computational Materials Science 43 (2008) 1207–1215

Finnis–Sinclair potentials. In fact the Finnis–Sinclair potentials also base on the second-moment approximation of tight-binding, however, their pair and electron density terms are expressed by polynomials [18]. To date, several potentials such as the Ackland and Thetford potential [19], Sutton–Chen [20] and extended F–S potential [21] can be considered as the modifications or extensions to the F–S potential. In addition, to improve the performance of empirical potentials in describing the interactions of atoms, Baskes et al. also made some modifications to these interatomic potentials and proposed modified EAM (MEAM) potentials, in which the angular factors are taken into account [22–24]. Whenas Vitek et al. proposed that the embedded energy could be expressed by two terms and both of them are functions of the electron density [25–27]. It has been shown by rule of thumb that the EAM and F–S potentials are applicable to studying the scientific issues related to fcc and bcc metals, whereas the SMA-TB potential has been extensively applied to fcc and hcp metals [28]. Hereby, a problem arises, i.e. how to describe the interatomic interactions in an hcp– bcc binary metal system. Concerning this problem, some researchers have attempted to construct EAM or SMA-TB potentials for hcp–bcc metal systems [28–31]. Nevertheless, these constructed potentials were not adequately validated and brought some problems [28]. One problem is related to the cutoff radius [32–34]. The cutoff radius defines the maximum distance of the interatomic interaction, beyond which the potential is truncated. In this case, proper treatment at the cutoff radius, i.e. incorporating truncation function is needed. Otherwise, the energy and force derived by these potential make jumps at the cutoff radius, and a large number of these events can spoil the energy conservation and lead to unphysical behaviors in simulations. Guellil and Adams have ever proposed a polynomial as the truncation function for EAM potentials (see Fig. 1) [32,35]. However, one sees that the cutoff problem has not yet been satisfyingly solved. Furthermore, the cutoff problem also exists in previously constructed potentials for some fcc– bcc binary metals, such as F–S and EAM potentials. In Fig. 1, the energy and force described by F–S potential have also been plotted, and one can clearly see that the jump at cutoff happens in F–S potential. In order to solve the cutoff-related problem, we recently proposed two empirical potentials for closed-packed transition metals and bcc transition metals, respectively [32,33]. Based on the newly proposed potential, we attempted to construct more relevant potentials for some hcp, fcc, bcc transition metals such as Zr, Hf, Cu, V, Nb and Ta and their binary alloys in the present study.

2. n-Body potential for transition metals 2.1. Construction of n-body potentials for Zr, Hf, Cu, V, Nb and Ta metals According to the formulism of SMA-TB, the energy of the d band is proportional to the square root of the second-moment of the density of states [4,16]. The total potential energy of an atom i can be written as Ei ¼

1X /ðrij Þ þ Fðqi Þ; 2 j6¼i

ð1Þ

in which pffiffiffiffi Fðqi Þ ¼  qi ;

ð2Þ

and qi ¼

X

wðr ij Þ:

ð3Þ

j6¼i

Here, rij is the distance between atoms i and j of the system at equilibrium state. For the sake of convenience, we ignore the precise physical argument here and may name / and w the pair function and density function, respectively. They can be written as 8 h  i > < A1 exp p1 rr0  1 ; 0  r  r cml ; /ðrÞ ¼ ð4aÞ  m h  i > : A2 rc1  r exp p2 rr0  1 ; rcml < r  r c1 ; r0 r0 and

wðrÞ ¼

h  i 8 r > > < B1 exp q1 r0  1 ;

0 < r  r cm2 :

 n h i > > : B2 rc2  r exp q ð r  1Þ ; 2 r0 r0 r0

rcm2 < r  r c2 ;

ð4bÞ

respectively. Here r0 is the distance between the nearest atoms. rc1 and rc2 are the cutoff radii of pair function and density function, respectively. rcm1 and rcm2 are the knot points of the segmented functions. m and n are two potential parameters and frequently are larger than 3, which can be designated by experience. Clearly, pair and density items and their high derivatives can continuously and smoothly go to zero at the cutoff radii rci, and rc2, respectively. A1, A2, B1, B2, p1, p2, q1 and q2, are potential parameters which should be determined. In practice, rc1, rc2 rcm1 and rcm2 can also be regarded as four potential parameters. One can see that there are total 14 parameters, denoted as {A1}, to describe the interaction among the fcc or hcp transition metals. The potential parameters can be determined through fitting the physical properties, such as lattice constants, cohesive energies and elastic constants to experimental values.

Table 1 Potential parameters of Zr, Hf and Cu metals

Fig. 1. Total energy and first derivatives as a function of lattice constants calculated from the F–S potential and EAM potential for bcc Ta (from Dai et al. [32]).

A1 (eV) A2 (eV) B1 (eV2) B2 (eV2) rcm1 (Å) rc1 (Å) rcm2 (Å) r2 (Å) r0 (Å) p1 p2 q1 q2 m n

Zr

Hf

Cu

0.3237 2.3721 4.3904 0.2571 3.5858 5.0355 4.0033 8.6898 3.1723 8.8624 0.1093 3.3925 0.0080 4 5

0.2465 1.0424 4.1993 0.0706 3.9582 5.3822 3.7030 9.7511 3.4289 11.3959 2.1761 6.1957 0.17311 4 6

0.1527 1.9297 1.4306 0.2289 2.4404 3.9497 3.3413 6.2061 2.5853 11.6683 4.8164 4.5780 0.0658 4 5

1209

J.H. Li et al. / Computational Materials Science 43 (2008) 1207–1215 Table 2 Lattice constants a and c (Å), cohesive energies EC (Å), bulk moduli, B0 and elastic constants Cij (Mbar) of hcp Zr, Hf and fcc Cu

a c EC B0 C11 C12 C13 C33 C55

Zr

Hf

Cu

3.18 3.23 5.19 5.15 6.25 6.19 0.89 0.83 1.42 1.43 0.67 0.73 0.55 0.65 1.54 1.64 0.26 0.32

3.29 3.19 5.38 5.05 6.43 6.41 1.05 1.09 1.77 1.88 0.76 0.77 0.61 0.66 1.93 1.69 0.36 0.56

3.66 3.62

3.49 3.48 1.37 1.37 1.69 1.68 1.21 1.22

hence eliminating the unphysical behaviors that might occur in molecular dynamics simulations. Moreover, the corresponding Rose et al. equations [37] are also plotted in Fig. 2. It is clear that the potential energies and the derivatives accord quite well with the Rose equations and almost overlapped in a considerably large range near the equilibrium point, indicating that the constructed potentials are relevant in describing the interatomic interaction among fcc and hcp transition metals. For bcc transition metals, such as V, Nb and Ta, the potential is also described by Eqs. (1)–(3). However, the pair and density functions are expressed by  m "  2  3  4 # r rc1 r r r r /ðrÞ ¼ A  1 þ c1 þ c2 þ c3 þ c4 ; r0 r0 r0 r0 r0 r0 0 < r < r c1 ;  n    r r c2 r wðrÞ ¼ B  exp b 1 ; r0 r  0 r0

0.79 0.76

The first lines are fitted and the second lines are measured in experiment (second line [36]).

The fitted potential parameters of hcp Zr, Hf and fcc Cu, are listed in Table 1. Table 2 shows the comparison between the fitted physical properties and that obtained in experiments [36]. It can be seen that these physical properties reproduced from the constructed potentials match well with the experimental results. To verify the validity of the constructed potential, Fig. 2 shows the energies and their derivatives of hcp Zr, Hf and fcc Cu derived from the constructed potential. One can see that both the energies and derivatives go smoothly to zero at the cutoff radii. Consequently, the ‘‘impulsive” contribution to the pressure are removed, and

ð5aÞ 0 < r < r c2 :

ð5bÞ

Here A, B, m, n, ci, b rci are potential parameters that can be determined by fitting to physical properties. Tables 3 and 4 show the fitted results of bcc Nb, Hf and Ta metals. Fig. 3 shows the energies and their derivatives of bcc Nb, Hf and Ta metals derived from the constructed potential. To further test the relevance of Eqs. (1)–(3) and Eq. (5), the surface energies of bcc Nb, V and Ta metals are calculated and compared both with experimental results and other method predictions. Foiles’ calculations show that the difference between the relaxed and unrelaxed surface energies is minor [38]. Therefore, the unrelaxed surface energies of three low-index, i.e. (1 1 0), (1 0 0) and (1 1 1) surfaces of bcc V, Nb and Ta are calculated by proposed potential and the calculated results are listed in Table 5. The surface energies calculated from other potentials and the experimental results [35,39,40] are also listed in Table 5. It can be seen Table 3 Potential parameters of V, Nb and Ta metals

8

Zr

4 0 -4 -8

E, dE/dx (eV)

8

Hf

4

A (eV) B (eV2) r0 (Å) rc1 (Å) rc2 (Å) c1 c2 c3 c4 b m n

V

Nb

Ta

350.644 0.2708 2.6241 4.6500 6.7000 3.7617 5.3347 3.3576 0.7870 1.8017 4 6

567.227 1.5940 2.8579 5.0200 6.7000 3.7319 5.2463 3.2679 0.7566 1.6823 4 6

441.485 1.7178 2.8579 5.0790 6.7000 3.7774 5.3924 3.4262 0.8145 0.1175 4 6

0 -4

Table 4 Lattice constants a (Å), cohesive energies EC (eV), bulk moduli B0 and elastic constants Cij (Mbar) of bcc V, Nb and Ta metals

-8 8

Cu

a

4 EC

0

E dE Rose

-4 -8 0.5

1.0

1.5

2.0

2.5

3.0

x = r/r0 Fig. 2. Potential energies E, their derivates dE as functions of x = r/r0 for the hcp Zr, Hf and fcc Cu derived from the constructed n-body potential together with the corresponding Rose equations.

B0 C11 C12 C55

V

Nb

Ta

3.03 3.03 5.31 5.31 1.57 1.57 2.29 2.29 1.21 1.21 0.44 0.44

3.30 3.30 7.57 7.57 1.72 1.71 2.47 2.47 1.35 1.35 0.29 0.29

3.30 3.30 8.10 8.10 1.94 1.94 2.66 2.66 1.58 1.58 0.87 0.87

The first lines are fitted and the second lines are measured in experiment (second line [36]).

1210

J.H. Li et al. / Computational Materials Science 43 (2008) 1207–1215 Table 5 Surface energies (J/m2) of bcc V, Nb and Ta metals

18.75 12.50

V

low-index surfaces

Methods

V

Nb

Ta

(1 1 0)

Present work EAM [35] MEAM [39] Experiment [40] Present work EAM [35] MEAM [39] Experiment [40] Present work MEAM [39]

1.676 1.683 2.636 1.827 1.929 1.831 2.778 2.596 2.161 2.931

1.792 1.807 2.490 1.859 2.101 1.968 2.715 2.628 2.343 2.923

1.909 1.800 2.778 2.019 2.223 1.990 3.035 2.868 2.480 3.247

6.25 0.00 (1 0 0)

E dE Rose

-6.25 -12.50

(1 1 1)

E, dE/dx (eV)

18.75

Nb

12.50

Table 6 Fitted cross potential of hcp–bcc Zr–Nb, Hf–Ta and fcc–bcc Cu–Nb systems

6.25 0.00

A1 (eV) A2 (eV) B1 (eV2) B2 (eV2) rcm1 (Å) rc1 (Å) rcm2 (Å) rc2 (Å) r0 (Å) p1 p2 q1 q2 m n

-6.25

18.75

Ta

12.50 6.25 0.00 -6.25

Zr–Nb

Hf–Ta

Cu–Nb

1.7577 9.6290 21.6441 2.0452 3.2181 5.1902 4.6120 7.3045 3.1396 6.3681 5.984105 7.0222 0.0259 4 6

0.8557 0.8448 10.1130 0.2596 0.1253 5.8061 5.2773 8.6407 3.0942 7.3910 0.0274 6.1956 1.5958 4 5

0.6831 3.0784 5.6772 4.5286 2.9275 4.7269 3.5786 5.5694 2.8056 8.5359 2.2991 7.0541 0.0076 5 5

-12.50 0.5

1.0

1.5

2.0

2.5

3.0

x Fig. 3. Potential energies E and their derivates dE as functions of x = r/r0 for the bcc V, Nb and Ta derived from long-range the long-range empirical potential and the corresponding Rose equations.

that for bcc V, Nb and Ta metals, the energy sequence of these faces is completely consistent with the predictions obtained by other potentials, i.e. the close-packed (1 1 0) face has the lowest energy, followed by the (1 0 0) and then the (1 1 1) faces, which are in good agreement with the experimental observations [41,42]. Moreover, the surface energies derived from present potentials are closer to experimental values than those obtained from the EAM potential, suggesting an apparent improvement comparing with those calculated from the empirical short-range potentials proposed previously. 2.2. Construction of n-body potentials for hcp–bcc and fcc–bcc binary metal systems In this section, we demonstrate the constructions of n-body potentials for hcp–bcc and fcc–bcc binary metal systems. For a mono-metal system, e.g. Zr metal, it is needed only one set potential parameter, {A1}Zr, to describe the interatomic interaction. However, for a binary metal system, e.g. Zr–Nb system, there should be three set potential parameters, e.g. {A1}Zr, {A1}Nb and {A1}Zr–Nb to characterize to describe the interaction between Zr–Zr, Nb–Nb and Zr–Nb atoms, respectively. Here {A1}Zr–Nb is often called as cross potential, which describes the interaction between atoms of different type. Generally, the cross potential {A1}Zr–Nb can also be determined by fitting some physical properties of Zr–Nb compounds. The physical properties of the compounds can either be measured in experiments or be acquired by ab initio calculations. In the present study, we intend to develop n-body potentials of the hcp–bcc structured Zr–Nb, Hf–Ta systems and fcc–bcc structured Cu–Nb system. The

four binary transition metal systems are all immiscible at equilibrium state, i.e. there is not any equilibrium compounds and hence there is not any physical properties can be used to fit the cross potentials. Consequently, the ab initio calculations are employed to assist constructing the interatomic n-body potential. The lattice constants, cohesive energies, elastic constants and bulk moduli of 6 L12 structured compounds [43,44], i.e. L12 Zr3Nb, ZrNb3, Hf3Ta, HfTa3, Cu3Nb, CuNb3 and 3 B2 structured compounds, i.e. B2 ZrNb, HfTa and CuNb, are first calculated by CASTEP [45,46], and then used in fitting procedure.

Table 7 Lattice constants a (Å), cohesive energies EC (eV), bulk moduli B0 and elastic constants Cij (Mbar) reproduced from the constructed potentials (first line) and acquired by ab initio calculations (second line) of L12 Zr3Nb, ZrNb3, Hf3Ta, HfTa3, Cu3Nb and CuNb3, and B2 ZrNb, HfTa and CuNb Structures

a

EC

B0

C11

C12

C44

L12 Zr3Nb

4.51 4.44 4.55 4.40 3.87 3.81

6.45 6.43 6.30 6.66 4.18 4.18

1.03 1.34 1.15 1.27 1.59 1.49

1.10 1.12 1.40 1.10 1.84 1.54

1.00 1.45 1.03 1.35 1.46 1.46

0.06 0.38 0.49 0.81 0.57 0.48

3.57 3.43 3.47 3.40 3.26 3.12 4.44 4.29 4.32 4.23 4.04 4.05

6.77 6.79 6.65 7.00 4.87 5.35 6.92 7.52 6.98 7.47 5.94 6.40

1.38 1.17 1.63 1.66 1.55 1.69 1.09 1.35 1.07 1.83 1.36 1.60

1.78 1.19 2.12 1.65 1.69 1.66 0.84 1.15 0.79 0.72 0.83 0.94

1.18 1.15 1.38 1.65 1.48 1.70 1.21 1.45 1.22 2.39 1.63 1.92

0.16 0.15 0.52 0.39 0.45 0.48 0.05 0.24 0.38 0.52 0.49 0.44

L12 Hf3Ta L12 Cu3Nb B2 ZrNb B2 HfTa B2 CuNb L12 ZrNb3 L12 HfTa3 L12 CuNb3

1211

J.H. Li et al. / Computational Materials Science 43 (2008) 1207–1215

parameters of cross potentials for Zr–Nb, Hf–Ta and Cu–Nb. Table 6 shows the fitted lattice constants, cohesive energies, bulk moduli and elastic constants of the metastable compounds. For comparison, the corresponding physical properties obtained by ab initio calculations are also listed in Table 7. One sees that the fitted physical properties match well with ab initio calculations. Fig. 4 shows the potential energies and their derivates as functions of the scale variable x = r/r0 for the L12 Zr3Nb, HfTa3 and B2 CuNb compounds. Fig. 4 shows that, for the potential energies well as their derivates, there is indeed no discontinuity in whole calculated range. It can also be seen that the result derived from the present potentials can keep smooth and matches well with that of Rose equation in whole calculated range, indicating that the present potential can reasonably describe the interatomic interaction of a system even far away from equilibrium state. To further test the relevance of the newly proposed potentials, some more physical properties of Zr, Hf, Cu, V, Nb, and Ta metals are calculated in the next Section.

18.75

L12 Zr3Nb

12.50 6.25 0.00

E dE Rose

-6.25 -12.50

E, dE/dx (ev)

12.50

L12 HfTa3

6.25

0.00

3. Physical properties of Zr, Hf, Cu, V, Nb and Ta metals derived from the proposed potentials

-6.25

3.1. Metastable structures of the Zr, Hf, Cu, V, Nb and Ta metals 12.50

The metastable structures of Zr, Hf, Cu, V, Nb and Ta metals are optimized and the corresponding lattice constants, cohesive energies, bulk moduli and elastic constants are then calculated from the potentials. Table 8 shows the calculated results. For comparison, the properties obtained by ab initio calculation are also listed in Table 8. It is noted that the lattice constants of these metastable structures are given by Wang et al. and the elastic constants as well as the bulk moduli are calculated by the present author with CASTEP employing PAW potentials [47]. Clearly, Table 8 shows that the reproduced properties especially the lattice constants, matching well with the ab initio calculation results. From the cohesive energies of these metastable structures, one can also see that the constructed potential can well reproduce the stable sequences of these transition metals. For example, the cohesive energies of the hcp and fcc V are 5.17 and 5.16 eV, respectively, both being lower than the cohesive energy 5.31 eV of bcc V. Based on the above calculations, it therefore can be concluded that the constructed potentials

B2 CuNb 6.25 0.00 -6.25 -12.50 0.5

1.0

1.5

2.0

2.5

3.0

x=r/r0 Fig. 4. Potential energies E and their derivates dE as functions x = r/r0for the L12 Zr3Nb, HfTa3 and B2 CuNb derived from the long-range empirical potential and the corresponding Rose equations.

Now we present the fitted results for the hcp–bcc Zr–Nb, Hf–Ta and fcc–bcc Cu–Nb binary metal systems. Table 5 shows the fitted

Table 8 Lattice constants a and c (Å), cohesive energies, EC (eV), bulk moduli and elastic constants, B0, and Cij (Mbar), of the metastable structures of V, Nb, Ta, Cu, Zr and Hf metals, reproduced from the constructed potentials (first line) comparing with that acquired by ab initio calculations (second line) hcp structure

fcc structure

bcc structure

V

Nb

Ta

Cu

V

Nb

Ta

Zr

Hf

Zr

Hf

Cu

2.61

2.84

2.86

2.59

3.70

4.04

4.06

4.49

4.65

3.57

3.70

2.91

2.88 4.73 5.24 7.36

2.90 4.77 5.16 7.91

2.56 4.22 4.20 3.48

3.81

4.23

4.22

4.53

4.47

3.56

3.54

2.86

EC

2.61 4.34 4.70 5.17

5.16

7.35

7.90

6.19

6.41

6.15

6.35

3.45

B0

– 1.23

– 1.26

– 1.50

– 1.37

– 1.31

– 1.28

– 1.34

– 0.88

– 1.05

– 0.84

– 1.01

– 1.32

C11

– 1.97

1.26 1.72

– 2.14

1.46 2.24

– 0.69

1.50 0.52

0.71

0.94 1.07

1.16 1.31

0.86 0.83

1.04 0.95

1.44 1.27

C12

– 1.23

– 1.37

– 1.22

2.62 1.02

1.62

0.50 1.67

1.66

1.30 0.78

1.86 0.91

0.92 0.85

0.24 1.04

0.91 1.35

– 0.67 – 1.94 – 0.18

– 0.85 0.46

1.03 0.84 0.77 2.42 2.75 0.43

–

2.50

0.76

0.81

0.84

1.45

1.71

0.33

– 0.60 – 2.57 – 0.15

0.76

0.49

0.90

0.50

0.66

0.57

0.79

0.94

–

0.43

–

0.31

–

0.49

–

0.65

0.65

0.20

0.36

0.90

a* c*

C13 C33 C55

*

The lattice constants are adopted from Ref. [47].

1212

J.H. Li et al. / Computational Materials Science 43 (2008) 1207–1215

are capable of in describing the interactions among the closedpacked transition metals and body-centered cubic transition metals. 3.2. Vacancy formation energies It is known that the interatomic interaction is environmentaldependent, i.e. the strength of the ‘‘individual bond” decreases as the local environment becomes very crowded. In order to examine the relevance of an n-body potential in describing the environmental dependence, the vacancy formation energy should be calculated and compared with experimental results. In the present study, an approximate method, proposed by Oh and Johnson [29], is employed to calculate the unrelaxed vacancy formation energy for EVunrelax 

N 1X /  qi F 0 ðqi Þ; 2 j6¼i ij

ð6Þ

where F0 (q) is the first derivative with respect to q. Table 9 gives the unrelaxed vacancy formation energies of bcc Nb, V, Ta, hcp Hf, Zr and fcc Cu metals calculated from the constructed potentials. For comparison, the experimental values and the results calculated by other researchers [48,49] are also listed in Table 9. One can see clearly that the vacancy formation energies calculated from the potentials agree well with the experimental observations at the quantity of about EC/3. 3.3. Melting points To estimate the melting points from the proposed potential, a series of molecular dynamics (MD) simulations at various temper-

Table 9 Unrelaxed vacancy formation energies calculated from the proposed potentials (first line) and vacancy formation energies obtained in experiments (second line)

Present work Experimental [48,49] Predicted by others researchers [48,49]

V

Nb

Ta

Zr

Hf

Cu

2.41

2.64 2.75 3.20 2.48 1.78 2.92

2.80 3.10 3.30 2.87 2.73 3.49

2.13 1.75 1.93 1.36

2.48 2.15 2.02

1.29 1.30 1.30

4000

Ta 3500

Nb

TM,MD (K)

3000

Hf 2500

V

Zr 2000

Cu T M,MD = 1.1T M,exp

1500

1500

atures are carried out for hcp Zr and Hf, fcc Cu, bcc V, Nb and Ta metals. The solid to liquid transitions, i.e. melting points, can be identified from the jumps in the temperature dependence of the total and potential energies obtained by MD simulations. Fig. 5 shows the comparison of the estimated melting points of the six studied metals in the present study and those obtained in experiments. It can be seen from the figure that the predicted melting points are averagely over-estimated by about 10%. One notes that in the MD simulations with periodic boundary conditions there are no possibilities for heterogeneous nucleation of the liquid phase and the crystal can be considerably overheated above the real melting points of the metals. Therefore, it can be concluded that the melting points determined by MD simulations in the present study are in well consistence with the experimental results.

2000

2500

3000

3500

T M,exp (K) Fig. 5. Comparison the melting points predicted by molecular dynamics simulation and experimental observation.

4. Formation and structure of Cu–Nb and Zr–Nb metallic glasses In the field of metallic glasses, i.e. amorphous alloys, one important and fundamental issue is to predict the composition range of a binary metal system, in which the metal glasses can be formed. Such composition range is often named glass-formation range of a metal system, which reflects the glass-formation ability of the system. Liu et al. proposed that the glass-formation range could be directly evaluated from interatomic potentials through MD simulation [50]. The principle of the method is to compare the relative stability of the supper-saturated solid solution versus its amorphous counterpart as a function of solute concentration, thus determining the two critical solid solubilities, beyond which the solid solutions become unstable and turn into the corresponding amorphous states. Accordingly, the glass-formation range of a system is thus deduced to be the composition range bounded by the two determined critical solid solubilities, in which an amorphous state is favored energetically comparing with its solid solution counterpart [50]. Based on the constructed potentials in the present study, MD simulations reveal that Cu–Nb and Zr–Nb metallic glasses could be formed within the composition range of 15–72 and 8–80 at.% Nb, respectively, matching well with experimental observations [51,52]. In the present study, we focus our intention on the dominating atomic packings of the Cu–Nb metallic glasses. In simulations, fcc, bcc and hcp structured solid solution models are constructed. For fcc and bcc models, the [1 0 0], [0 1 0] and [0 0 1] crystal directions are parallel to the x-, y- and z-axes, respectively. For hcp models, the [1 0 0], [1 2 0] and [0 0 1] crystalline directions are parallel to the x-, y- and z-axes, respectively. There are total 2048, 2000 and 2352 atoms in the fcc, bcc and hcp simulation models, respectively. In setting simulation models, the solute atoms are added by random substituting the desired number of solvent atoms, obtaining the initial state of solid solutions with desired compositions. For the three solid solution models, the minimum image criterion must be satisfied and the periodic boundary conditions are adopted in x-, y- and z-directions. A series of MD simulations are then conducted for those models of different compositions at 300 K and zero pressure with Parrinello–Rahman constant pressure scheme [53]. The equations of motion are solved using the second-order fourth-values predictor-corrector algorithm of Gear with a time step 10  1015 s [54]. The simulation time is about 50,000 to 150,000 MD time steps for each run to reach a relatively equilibrium (metastable) state, at which all the related dynamic variables show no secular variation. If one wants to estimate the melting point of a solid solution, one should perform a series of MD simulations at different temperatures and obtain the energy–temperature and volume–temperature curves. The jumps in these curves indicate the solid to liquid transition, from which the melting point can thus be determined.

1213

J.H. Li et al. / Computational Materials Science 43 (2008) 1207–1215

ities obtained in the present study. It can be seen from Fig. 6a that the coordination polyhedra of CN = 12 dominates the atomic packings, corresponding to the bcc structures, in the models dissolving 20 and 25 at.% Cu. The fractions of the dominating atomic packings are both larger than 80%. Once the concentration exceeds the critical solid solubility, e.g. in the models dissolving 30 and 35 at.% Cu, the dominating coordination polyhedra are CN = 14, 15 and 13, exhibiting the typical characteristic of metallic glasses. For the solid solution obtained for fcc simulation models, one notes an interesting phenomenon in Fig. 6b, i.e. the dominating coordination polyhedra is CN = 14 even for pure Cu. In fact, in an ideal fcc structure, the Voronoi polyhedron is dodecahedral. However, a very small displacement of atoms in this structure may destroy the dodecahedron and form an icosidihedron, an icositetrahedron, an icosihexahedron or other type Voronoi polyhedrons, which indicates that the Voronoi method are not good for fcc structured alloys [56]. Even so, Fig. 6b still shows the slight differences between faced-centered cubic and disordered structures. For example, in the alloys obtained from the simulation models dissolving 0, 5 and 10 at.% Nb, the fractions of the dominating coordination polyhedra (CN = 12) are larger than 40%, which corresponds to fcc structures. When the concentration of Nb exceeds the critical solid solubility, the Cu–Nb alloys obtained in the models dissolving 15, 20 and 25 at.% Cu shows typical amorphous characteristics. The dominating coordination polyhedra are CN = 14, 15 and 13. Fig. 7 shows the correlation between the composition and fractions of coordination polyhedra in Cu–Nb alloys obtained by the present simulations. It can be seen that the Cu–Nb metallic glasses could

It is known that the coordination numbers and/or coordination polyhedrons are usually employed to characterize the atomic configuration of amorphous alloys. The coordination polyhedrons or coordination numbers can be defined by and derived from the Voronoi polyhedrons (i.e. generalized Wigner–Seitz cells) that can be precisely and uniquely determined by the Voronoi method [55]. Generally, the faces, vertices and edges of a Voronoi polyhedron in a disordered structure correspond to the vertices, faces and edges of a coordination polyhedron, respectively. Consequently, the structural change in the solid solution models is characterized using the Voronoi analyses. Table 10 shows the Voronoi polyhedra and corresponding coordination polyhedra observed in Cu–Nb metallic glasses in the present study. One should note that the number of the faces and edges of the coordination polyhedron derived from the correspondence is maximum possible values and may be reduced in some symmetrical structures. For example, there are 14 faces and 24 vertices in the Voronoi polyhedron of an ideal bcc structure, however, the coordination polyhedron is dodecahedral. The Voronoi polyhedron of an ideal fcc structure is dodecahedral while the numbers of faces, vertices and edges of the coordination polyhedron are 14, 12 and 24, respectively, but not 20, 12 and 30 given in Table 10. Nevertheless the number of vertices of the coordination polyhedra, i.e. the coordination number (CN), always equal to the number of the faces of the Voronoi polyhedron. Therefore, we would use the coordination numbers to denote the atomic packings of metallic glasses. Fig. 6 shows the fractions of coordination polyhedra in Cu–Nb alloys with compositions in the vicinity of the critical solid solubil-

Table 10 Voronoi polyhedra and corresponding coordination polyhedra observed in Cu–Nb metallic glasses in the present study Voronoi polyhedra

Enneahedron Decahedron Hendecahedron Dodecahedron Tridecahedron Tetradecahedron Pentadecahedron Hexadecahedron Heptadecahedron Octadecahedron Enneadecahedron

Coordination polyhedra

nf

nv

ne

9 10 11 12 13 14 15 16 17 18 19

14 16 18 20 22 24 26 28 30 32 34

21 24 27 30 33 36 39 42 45 48 51

Tetradecahedron Hexadecahedron Octadecahedron Icosahedron Icosidihedron Icositetrahedron Icosihexahedron Icosioctahedron Triacontahedron Triacontadihedron Triacontatetrahedron

Fractions

nf

nv

ne

(%)

14 16 18 20 22 24 26 28 30 32 34

9 10 11 12 13 14 15 16 17 18 19

21 24 27 30 33 36 39 42 45 48 51

<0.1 <0.1 0.5 7.7 20.7 33.3 26.2 9.5 1.9 <0.1 <0.1

nf, nv and ne indicate the numbers of faces, vertices and edges of the polyhedra. One notes that nf and ne of the coordination polyhedra may reduce in some special cases.

a

b

0.5

1.0 20% Cu 25% Cu 30% Cu 35% Cu

bcc

0.6

0.4

fcc 75% Cu 80% Cu 85% Cu 90% Cu 95% Cu 100% Cu

amor.

0.4

amor.

0.3

0.2

0.2

0.1

0.0

0.0 12

14

16

CN

18

20

12

14

16

18

Fractions (%)

Fractions (%)

0.8

20

CN

Fig. 6. Fractions of coordination polyhedra in Cu–Nb alloys with compositions in the vicinity of the critical solid solubilities: (a) obtained from bcc structured simulation models and (b) obtained from fcc structured simulation models.

1214

J.H. Li et al. / Computational Materials Science 43 (2008) 1207–1215

a

1.0

b

Fractions (%)

GFR

CN= 12 13 14 15 16 17

0.6 0.4

0.4

Fractions (%)

0.6 0.8

0.2

0.2 0.0

0.0 0

10

20

30

40

50 50

60

70

80

90

100

C (Cu%)

C (Cu%)

Fig. 7. Correlation between the composition and fractions of coordination polyhedra in Cu–Nb alloys obtained by MD simulations: (a) obtained from bcc structured simulation models and (b) obtained from fcc structured simulation models.

lated from the interatomic potential can go smoothly to zero at cutoff radius, apparently avoiding the unphysical behaviors, such as the abnormal fluctuation of energy and pressure maybe emerging in computer simulations. The constructed potentials have been applied to study the amorphization of the Cu–Nb and Zr–Nb systems by molecular dynamics simulation. It is found that Cu–Nb and Zr–Nb metallic glasses could be formed within the composition range of about 15–72 and 8–80 at.% Nb, respectively, matching well with experimental observations. Voronoi analyses reveal that the dominating atomic packings in Cu–Nb metallic glasses are the icositetrahedron (CN = 14), icosihexahedron (CN = 15) and icosidihedron (CN = 13) with fractions of 33%, 26% and 21%, respectively.

0.40 0.35

bcc model (30% Cu) fcc model (30% Cu) bcc model (50% Cu) fcc model (50% Cu) bcc model (70% Cu) fcc model (70% Cu)

Fractions (%)

0.30 0.25 0.20 0.15 0.10 0.05

Acknowledgments

0.00 6

8

10

12

14

16

18

20

CN Fig. 8. Fractions of coordination polyhedra in Cu–Nb metallic glasses with different compositions and obtained from different simulation models.

The authors are grateful to the financial support from the National Natural Science Foundation of China (50531040), The Ministry of Science and Technology of China (2006CB605201), and the Administration of Tsinghua University. References

be formed when the composition falls into the range of about 15– 72 at.% Nb. The glass-formation range determined in the present study matches well with the experiment results [51,57]. Fig. 7 also suggests that the fractions sequences of coordination polyhedra are almost identical for the Cu–Nb amorphous alloys with different compositions. To further confirm this point, the Cu–Nb amorphous alloys with 30, 50 70 at.% Cu obtained from fcc and bcc simulation models are analyzed by Voronoi method. The results are shown in Fig. 8, which suggests that the dominating coordination polyhedra in the six Cu–Nb metallic glasses are all CN = 14 (icositetrahedron), 15 (icosihexahedron) and 13 (icosidihedron) with fractions of 33%, 26% and 21%, respectively. 5. Conclusion In the framework of second-moment approximation of the tight-binding theory, the empirical n-body potentials are constructed for close-packed transition metals, i.e. hcp Zr, Hf, fcc Cu and bcc-structured V, Nb, Ta metals. With the assistance of ab initio calculation, the interatomic potentials are constructed for hcp–bcc Zr–Nb, Hf–Ta and fcc–bcc Cu–Nb binary metal systems. It turns out that, for the hcp, fcc and bcc transition metals and their alloys, the constructed potentials can correctly predict the vacancy formation energies and the stable sequence of different structures of these metals. More importantly, both their energies and pressures calcu-

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

J.E. Lennard-Jones, Proc. Royal Soc. 109 (1924) 584. Erkocß Sßakir, Phys. Rep. 278 (1997) 79. M.S. Daw, S.M. Foiles, M.I. Baskes, Mater. Sci. Rep. 10 (1993) 251. V. Rosato, M. Guillope, B. Legrand, Philos. Mag. A 59 (1989) 321. M. Finnis, Interatomic Force in Condensed Matter, Oxford, 2003. M.S. Daw, M.I. Baskes, Phys. Rev. Lett. 50 (1983) 1285. S.M. Foiles, M.I. Baskes, M.S. Daw, Phys. Rev. B 33 (1986) 7983. S.M. Foiles, M.S. Daw, Phys. Rev. B 38 (1988) 12643. M.S. Daw, Phys. Rev. B 39 (1989) 7441. R.A. Johnson, Phys. Rev. B 37 (1988) 3924. R.A. Johnson, D.J. Oh, J. Mater. Res. 4 (1989) 1195. J. Cai, Y.Y. Ye, Phys. Rev. B 54 (1996) 8398. H.N.G. Wadley, X. Zhou, R.A. Johnson, M. Neurock, Prog. Mater. Sci. 46 (2001) 329. Y. Mishin, M.J. Mehl, D.A. Papaconstantopoulos, Acta Mater. 53 (2005) 4029. P.L. Williams, Y. Mishin, J.C. Hamilton, Model. Simul. Mater. Sci. Eng. 14 (2006) 817. F. Willaime, C. Massobrio, Phys. Rev. B 43 (1991) 11653. F. Cleri, V. Rosato, Phys. Rev. B 48 (1993) 22. M.W. Finnis, J.E. Sinclair, Philos. Mag. A 50 (1984) 45. G.J. Ackland, R. Thetford, Philos. Mag. A 56 (1987) 15. A.P. Sutton, J. Chen, Philos. Mag. Lett. 61 (1990) 139. X.D. Dai, Y. Kong, J.H. Li, B.X. Liu, J. Phys. Condens. Matter 18 (2006) 4527. M.I. Baskes, Phys. Rev. Lett. 59 (1987) 2666. Y.F. Ouyang, B.W. Zhang, Z.P. Jin, S.Z. Liao, Z. Phys. B 87 (1996) 802. Y.M. Kim, B.J. Lee, M.I. Baskes, Phys. Rev. B 74 (2006) 014101. M. Igarashi, M. Khantha, V. Vitek, Philos. Mag. B 63 (1991) 603. W.Y. Hu, X.L. Shu, B.W. Zhang, Comput. Mater. Sci. 23 (2002) 175. Y.R. Wu, W.Y. Hu, S.S. Han, J. Alloy. Compd. 420 (2006) 83. R.F. Zhang, Y. Kong, B.X. Liu, Phys. Rev. B 71 (2005) 214102. D.J. Oh, R.A. Johnson, J. Mater. Res. 3 (1988) 471. R. Pasianot, E.J. Savino, Phys. Rev. B 45 (1992) 12704.

J.H. Li et al. / Computational Materials Science 43 (2008) 1207–1215 [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]

D. Nguyen-Manh, S.L. Dudarev, Mater. Sci. Eng. A 423 (2006) 74. X.D. Dai, J.H. Li, Y. Kong, Phys. Rev. B 75 (2007) 052102. J.H. Li, X.D. Dai, T.L. Wang, B.X. Liu, J. Phys. Condens. Matter 19 (2007) 086228. D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic Press, New York, 2002. A.M. Guellil, J.B. Adams, J. Mater. Res. Res. 7 (1992) 639. David R. Lide, CRC Handbook of Chemistry and Physics: A ready-reference Book of Chemical and Physical Data, CRC Press, Boca Raton, FL, 2004. F. Guinea, J.H. Rose, J.R. Smith, J. Ferrante, Appl. Phys. Lett. 44 (1984) 53. S.M. Foiles, Surf. Sci. 191 (1987) 779. B.J. Lee, M.I. Baskes, H. Kim, Y.K. Cho, Phys. Rev. B 64 (2001) 184102. W.R. Tyson, W.A. Miller, Surf. Sci. 62 (1977) 267. B.E. Sundquist, Acta Metall. 12 (1964) 67. H.E. Grenga, R. Kumar, Surf. Sci. 61 (1976) 283. W. Pearson, A Handbook of Lattice Spacings and Structures of Metals and Alloys, Pergamon Press, London, 1958. C.J. Smithells, Smithells Metals Reference Book, Butterworth-Heinemann, Boston, 1992.

1215

[45] M.D. Segall, P.L.D. Lindan, M.J. Probert, C.J. Pickard, P.J. Hasnip, S.J. Clark, M.C. Payne, J. Phys. Condens. Matter. 14 (2002) 2717. [46] B.X. Liu, W.S. Lai, Z.J. Zhang, Adv. Phys. 50 (2001) 367. [47] Y. Wang, S. Curtarolo, C. Jiang, R. Arroyave, T. Wang, G. Ceder, L.Q. Chen, Z.K. Liu, Calphad-comput. Coup. Phase Diag. Thermochem. 28 (2004) 79. [48] Avinc Ahmet, Dimitrov Ventzislav Ivanov, Comput. Mater. Sci. 13 (1999) 211. [49] S.L. Chen, J.C. Xu, H.S. Zhang, Comput. Mater. Sci. 29 (2004) 428. [50] B.X. Liu, W.S. Lai, Q. Zhang, Mater. Sci. Eng. R 244 (2000) 1. [51] T.L. Wang, J.H. Li, K.P. Ta, B.X. Liu, Scripta Mater. 57 (2007) 157. [52] S.H. Liang, J.H. Li, B.X. Liu, Comput. Mater. Sci. 42 (2008) 550. [53] M. Parrinello, A. Rahman, J. Appl. Phys. 52 (1981) 7182. [54] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford, Clarendon, UK, 1981. [55] V.S. Stepanyuk, A. Szasz, A.A. Katsnelson, O.S. Trushin, H. Müller, H. Kirchmayr, J. Non-Cryst. Solids 159 (1993) 80. [56] Faken Daniel, Jonsson Hannes, Comput. Mater. Sci. 2 (1994) 279. [57] C. Michaelsen, C. Gente, R. Bormann, J. Appl. Phys. 81 (1997) 6024.