An analytic n-body potential for bcc Iron

NIM B Beam Interactions with Materials & Atoms

Nuclear Instruments and Methods in Physics Research B 255 (2007) 37–40 www.elsevier.com/locate/nimb

An analytic n-body potential for bcc Iron V. Pontikis a

a,*

, V. Russier b, J. Wallenius

c

Commissariat a` l’Energie Atomique, DRECAM/LSI, CE de Saclay, Building 524, Room 40B, 91191 Gif-sur-Yvette Cedex, France b Centre d’Etudes de Chimie Me´tallurgique, CNRS UPR2801, 94407 Vitry-sur-Seine, France c Royal Institute of Technology, Department of Nuclear and Reactor Physics, Stockholm, Sweden Available online 2 January 2007

Abstract We have developed an analytic n-body phenomenological potential for bcc iron made of two electron-density functionals representing repulsion via the Thomas–Fermi free-electron gas kinetic energy term and attraction via a square root functional similar to the second moment approximation of the tight-binding scheme. Electron-density is given by radial, hydrogen-like orbitals with effective charges taken as adjustable parameters fitted on experimental and ab-initio data. Although the set of adjustable parameters is small, prediction of static and dynamical properties of iron is in excellent agreement with the experiments. Advantages and shortcomings of this model are discussed with reference to published works. Ó 2006 Elsevier B.V. All rights reserved. PACS: 34.20.Cf; 61.50.Ah; 61.50.Lt; 65.40.De; 65.40.Gr Keywords: n-body potentials; Atomistic simulation; Static and dynamical properties

1. Introduction The simplified representation of atomic cohesion via phenomenological potentials is still an attractive methodology to predict equilibrium and non-equilibrium properties of metals and alloys via large-scale atomistic simulations. This is how macroscopic properties of materials are obtained from electronic structure, atomistic and mesoscopic calculations, each level taking input data from the next smaller scale [1]. Customarily, phenomenological potentials for metals and alloys are made of a physically motivated, n-body cohesive term and an empirical repulsive term schematically representing coulombic interactions and effects of Pauli’s exclusion principle. One serious obstacle to the transferability of such models is the empirical character of the repulsive term, which does not contribute satisfactorily to the energy and its derivatives whenever a single reference state is used in fitting the parameters of *

Corresponding author. Tel.: +33 1 6908 2904; fax: +33 1 6908 2199. E-mail address: [email protected] (V. Pontikis).

0168-583X/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.nimb.2006.11.050

the potential. The common practice adopted in the literature to bypass the deficiency of this term is expanding it in as much as needed functions so that to satisfy the requirements of several reference states, either experimental or resulting from ab-initio calculations. Thus the number of adjustable parameters rapidly increases without improving the physical content of the model [2,3]. This is a strong motivation to describe the repulsion by a phenomenological yet physically motivated term a strategy also expected to reduce the number of adjustable parameters. In the present work we have implemented a n-body phenomenological potential for iron such that: (i) the repulsive term is the Thomas–Fermi free-electron gas functional applying to the electronic density of 4s valence electrons (ii) the attractive term is represented by a square root functional of the electronic density of 3d valence electrons and (iii) electronic density in the latter case is calculated by using solely the radial part of the effective atomic wave function. We show that within this scheme, a reduced set of adjustable parameters leads to a remarkably good agreement between calculated and experimental static and

38

V. Pontikis et al. / Nucl. Instr. and Meth. in Phys. Res. B 255 (2007) 37–40

dynamical properties. In Section 2, the model and computational details are given whereas Section 3 and 4 are respectively devoted to the presentation of the results, a short discussion and some conclusive remarks. 2. Model and computations The total energy of N iron atoms is written as the sum of site energies of point particles: X E¼ U i; ð1aÞ i

Ui ¼ A

X j6¼i

!5=3 q4s ðrij Þ

n

X

!1=2 q3d ðrij Þ

;

ð1bÞ

j6¼i

with the electron-densities of 4s and 3d electrons given by:  2 ð2aÞ q4s ðrij Þ ¼ W4s ðrij Þ ;  2 q3d ðrij Þ ¼ W3d ðrij Þ ð2bÞ and hydrogen-like radial wave functions, Wa(r), written as follows:  2 1 24  18Z 4s r þ 3 Z 4s r 96   3 !   Z r Z r  4s ð3aÞ exp  4s ; 2 4   2   1 2Z 3d r Z 3d r W3d ðrÞ ¼ pffiffiffiffiffi exp  ð3bÞ 3 3 9 30 with A, n, Z 4s and Z 3d , taken as adjustable parameters and r the interatomic distances expressed in atomic units. However, the space decay of the above listed atomic wave functions is incompatible with the usually short-range interactions in transition metals. Therefore, the range of the electronic densities, q4s and q3d (Eq. 1b), is explicitly limited via a Fermi–Dirac step function, f(r, rc, e) = 1/(1+ exp (e(r/rc  1))), acting as a multiplicative factor. The resulting phenomenological potential is a simple analytic model of cohesion with six adjustable parameters to be fitted on experimental and/or ab-initio data. It is worth noting that repulsion in this model has the functional dependence of the kinetic energy in a Thomas–Fermi free electron gas and that accordingly the electronic density entering this term for iron is that of 4s electrons. Moreover, core electrons and gradient of the electronic density contributions to the cohesion are not considered. These become important when atoms approach each other at very short distances, a situation in which this model would fail describing the very short-range repulsion. On the other hand, cohesion in transition metals is mainly accounted for by d-electrons [4]. Therefore, the attractive term consists in a square-root functional of the electronic density, by analogy to the second moment approximation of tight-binding and 3d effective atomic wave functions enter in this term. The model parameters have been fitted to a selected set of experimental data by using MERLIN, a multidimenW4s ðrÞ ¼

sional minimization package [5]. With an optimal set of parameters, the potential was tested by computing phonon dispersion, the temperature dependence of the lattice constant and of the atomic mean square displacements (MSD) and the formation energies of interstitials. Energy minimizations have been performed via a quasi-dynamic procedure [6] whereas molecular dynamics (MD) in the canonical and isothermal–isobaric ensembles have been made according to the Nose´–Andersen scheme [7,8]. MSD’s have been computed over constant temperature equilibrium trajectories, which lasted 5000 time steps (dt = 1015 s). 3. Results Table 1 collects the numerical values of the potential parameters whereas in Table 2 experimental data used in the fit are compared with their calculated counterparts, ab-initio data [3] and data obtained with an embedded atom method (EAM) potential from the literature [2,3]. It can be seen from Table 2 that the agreement is excellent between calculated and experimental target values for all fundamental state properties at T = 0 K. Fig. 1 shows further that the agreement between model and experiment [9,10] is also good for the phonon dispersion, though small differences subsist at the limits of the Brillouin zones. This implies that the model should correctly predict the free energy change of a-iron with temperature. The agreement between calculations and experiments [11] is also satisfactory for the lattice constant and MSD’s as a function of the temperature (Figs. 2 and 3). This is not always true for other phenomenological models including those leading to acceptable values for the linear expansion coefficient. These predict lattice constant values significantly different from experimental ones with as a main consequence MSD’s far off their experimental counterparts. It is worth noticing that the a-iron potential recently developed by Wallenius et al. [3] makes exception to this trend at the extra cost of violating the experimental equation of state at short distances. The formation energies of point defects are another crucial test of transferability of phenomenological potentials. Concerning the vacancy formation energy, the present model behaves satisfactorily (Table 2). However, it seriously fails reproducing the formation energy of interstitials, which is found ranging from 0.6 eV (split h1 1 1i) to 1.7 eV (split h1 1 0i). This behavior is fully consistent with the remarks made above in Section 2: at short distances the electronic density deviates from the homogeneous disTable 1 Values of the potential parameters and of the spline coefficients, As, rs ˚ 3) rs A (eV) n (eV) Z Z e rc As (eV/A 1011

147.9

4s

3d

3.15

0.507

15.47

1.035

2137.33

0.8

Effective charges, Z 4s and Z 3d are in elementary charge units whereas distances, rc and rs, are in units of the lattice constant at T = 0 K, a0 = 0.286 nm.

V. Pontikis et al. / Nucl. Instr. and Meth. in Phys. Res. B 255 (2007) 37–40

39

Table 2 Properties of a-Fe calculated by using the potential developed in this work This work

Experiment

VASP

EAM [2,3]

a0 (nm) B (GPa) C 0 (GPa) C44 (GPa) Ecoh (eV) Evf (eV) Efcc  Ebcc (eV) Efh110i

0.286 172.7 52.8 121.8 4.289 1.78a 0.03 5.83a

0.286 [14,15] 173.1 [14] 52.5 [14] 121.8 [14] 4.28 [16] 2 ± 0.2 [17] 0.05 [18] 3.0–12.0 [19,20]

0.283 [3] – – – – 2.12 [3] 0.09 [3] 3.4–4.0 [3,12,13]

0.2855 178 49.2 116 4.01 1.73 0.12 3.45

Efh110i  Efh111i

0.11

–

0.7, 0.67 [3,12,13]

0.5

The comparison is made with ab-initio calculations [3,12,13], data obtained by using an EAM potential from the literature [2,3] and experimental values extrapolated at T = 0 K [14–18]. a Relaxed value.

Γ [ζ00]

[ζζζ] Γ [ζζ0 ]

H

10

L

L

2.9

T

6

T

T

4

T

α-Fe lattice constant (Å)

L

8

ω (THz)

2.91

2

1

2.89

2.88

2.87

2.86

2

2.85

0 0

0.5

1

1.5

2

0

2.5

400

1200

T (K)

k (2π/a) Fig. 1. Phonon dispersion: experiment (open squares, triangles and circles, [9,10]) and lattice dynamics (full lines).

tribution and the absence of core electrons leads to the observed failure. Owing to the good quality offered by this model for the properties of the fundamental state, the classical solution found in the literature has been adapted to the present case i.e. a repulsive spline acting at short distances has been added to the potential: U ðrÞ ¼ As ðrs  rÞ3 H ðrs  rÞ

800

ð4Þ

with a cutoff radius, rs = 0.8a0, entering the definition of Heaviside’s step function, H(rs  r), significantly smaller than the distance between first neighbors, d1NN  0.866a0. It is worth mentioning that this additional term does not affect fundamental state properties over the entire domain of atomic densities of solid and liquid iron at low pressures. An optimal choice of the spline coefficients is given in Table 1 leading to formation energies of split interstitials qualitatively consistent with experiment (Table 2) though

Fig. 2. Lattice constant as a function of the temperature: experiment (full circles and triangles, [14,15]) and MD computations (full diamonds, the line is just a guide for the eye).

their magnitudes and relative stability are in contrasting difference with electronic structure calculations [3,12,13]. 4. Discussion and perspectives Among early attempts to model cohesion in transition metals, the tight-binding approach at the second moment approximation has revealed successful and remarkably simple for describing trends of the energy and its derivatives in compact lattices as a function of the number of d-electrons [21,22]. The case of body centered cubic transition (bcc) metals is similar, provided at least second neighbors are explicitly considered [4,22]. The present model is a straightforward extension to bcc metals of the second moment approximation of tight-binding augmented by an n-body repulsive interaction related to the kinetic energy of s-electrons. Its main quality consists in a remarkably

40

V. Pontikis et al. / Nucl. Instr. and Meth. in Phys. Res. B 255 (2007) 37–40

tronic density we have here adopted can also closely fit the electronic density as obtained from ab-initio calculations and leads thereby to better n-body potentials than one can obtain via the widespread empirical procedures. Current work focuses on improving the model and applying it to bcc chromium and the iron–chromium alloy.

0.01

-16

2

(x 10 cm )

0.02

x

2

References

0.01

0.00

0

200

400

600

800

1000

T (K) Fig. 3. Atomic mean square displacements as a function of the temperature: experiment (open squares, the line is just a guide for the eye) and MD computations (full circles).

good description of ground state properties of perfect bcc metals though simple, analytic and still with a small number of adjustable parameters by comparison with other models found in the literature [2]. It is worth noting that the present model does not account satisfactorily for the formation energies and the stability of split interstitials. This is the justification for adding the spline term, a wellknown remedy removing this limitation [2,3]. However, no attempt has been made here of further refining this additional repulsion term given that such an effort would not increase the physical content of the model. Indeed in iron, the stability order of interstitials and absolute values of formation energies are suggested to correlate with magnetism and thus can only indirectly be reflected in the values of parameters of empirical terms [23,13]. A last comment is of interest about the values taken in this model by the effective charges Z4s and Z3d: though Z4s is close to the value 3.75 prescribed by Slater’s rules [24], Z3d is far off this prescription, 6.25, emphasizing on that no physical meaning should be given to the values here obtained. Nevertheless, in a similar in spirit approach Mitev et al. [25] have recently shown that the representation scheme of the elec-

[1] H.O. Kirchner, L.P. Kubin, V. Pontikis (Eds.), Computer Simulation in Materials Science: Nano-/Meso-/Macroscopic Space and Time Scales, NATO ASI Series E, Vol. 308, 1995. [2] M.I. Mendelev, S. Han, D.J. Srolovitz, G.J. Ackland, D.Y. Sun, M. Asta, Philos. Mag. 83 (2003) 3977. [3] J. Wallenius, P. Olsson, C. Lagerstedt, Nucl. Instr. and Meth. B 228 (2005) 122. [4] F. Ducastelle, J. Phys. 31 (1970) 1055. [5] G.A. Evangelakis, J.P. Rizos, I.E. Lagaris, I.N. Demetropoulos, Comput. Phys. Commun. 46 (1987) 401. [6] J.R. Beeler, G.L. Kulcinski, Interatomic Potentials and Simulation of Lattice Defects, in: P.C. Gehlen, J.R. Beeler, R.I. Jaffee (Eds.), Plenum Press, New York, 1972, p. 735. [7] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [8] S. Nose´, Mol. Phys. 52 (1984) 255. [9] B.N. Brockhouse, H.E. Abou-Helal, E.D. Hallman, Solid State Commun. 5 (1967) 221. [10] V.J. Minkievicz, G. Shirane, R. Nathans, Phys. Rev. 162 (1967) 528. [11] E.A. Owen, E.W. Evans, Br. J. Appl. Phys. 18 (1967) 611. [12] C.C. Fu, F. Willaime, P. Ordejon, Phys. Rev. Lett. 92 (2004) 175503. [13] C. Domain, C. Becquart, Phys. Rev. B 65 (2001) 024103. [14] G. Simmons, H. Wang, Single Crystal Elastic Constants and Calculated Aggregate Properties: A Handbook, The MIT Press, Cambridge, 1971. [15] Z.S. Basinski, W. Hume-Rothery, A.L. Sutton, Proc. Royal Soc. A 229 (1955) 459. [16] C. Kittel, Introduction to Solid State Physics, fifth ed., Wiley, 1976, Table 3. [17] L.D. Schepper, D. Segers, L. Dorikens-Vanpraet, M. Dorikens, G. Knuyt, L. Stals, P. Moser, Phys. Rev. B 27 (1983) 5257. [18] W. Bendick, W. Pepperhof, Acta Metall. 30 (1982) 679. [19] P. Moser, Mem. Sci. Rev. Metall. 63 (1966) 431. [20] H. Bilger, V. Hivert, J. Verdone, J. Leveque, J. Soulie, in: International Conference on vacancies and interstitials in metals, Kernforshungsanlage Juelich, 1968, p. 751. [21] V. Rosato, M. Guillope´, B. Legrand, Philos. Mag. A 59 (1989) 321. [22] F. Ducastelle, J. Phys. C 7 (1974) 79. [23] D. Nguyen-Manh, A.P. Horsfield, S.L. Dudarev, Phys. Rev. B 73 (2006) 020101. [24] J.C. Slater, Phys. Rev. 36 (1930) 57. [25] P. Mitev, G.A. Evangelakis, E. Kaxiras, Model. Simul. Mater. Sci. Eng. 14 (2006) 721.