Nuclear Instruments and Methods in Physics Research B 228 (2005) 245–249 www.elsevier.com/locate/nimb
Coupling of MC and MD techniques for the calculation of vacancy cluster binding energies D. Kulikov a
a,b
, L. Malerba
c,* ,
M. Hou
a
Physique des Solides Irradie´s et des nanostructures CP234, Universite´ Libre de Bruxelles, Bd du Triomphe, B-1050 Brussels, Belgium b A.F. Ioffe Physico-Technical Institute of RAS, Polytechnicheskaya Str. 26, 194021 St. Petersburg, Russia c SCKCEN, Nuclear Research Centre, Reactor Materials Research Unit, Boeretang 200, B-2400 Mol, Belgium
Abstract A combination of Metropolis Monte Carlo (MMC) and molecular dynamics (MD) techniques has been used to systematically sample configurations of vacancy clusters (VC) in zirconium and calculate the corresponding binding energies as a function of cluster size. The application of this procedure allowed the interpolation of expressions for the binding energies, which are of immediate use in kinetic Monte Carlo (KMC) or rate theory (RT) models. The results of the application of the procedure are presented, analysed and discussed. 2004 Elsevier B.V. All rights reserved.
1. Introduction Kinetic Monte Carlo (KMC) techniques [1], as well as rate theory (RT) [2], are efficient tools for the study of the long-term evolution of radiation damage. The main problem for the application of these models is that the knowledge of a large number of fundamental magnitudes, such as defect migration energies, reaction distances and cluster binding energies, is a necessary prerequisite. The binding energies of vacancy clusters are particu*
Corresponding author. Tel.: +32 14 333090; fax: +32 14 321216. E-mail address:
[email protected] (L. Malerba).
larly important because they will decide the rate of dissociation of the clusters as a function of temperature, thereby governing their nucleation and growth rate. Unfortunately, these are not magnitudes that can be experimentally determined and must be calculated. Semi-empirical methods can be used for such purpose. Calculations of this type are heavy, because of the large number of a priori possible cluster configurations, even at constant size. In this work a computational procedure to systematically calculate the formation energies and corresponding binding energies is proposed and applied to the case study of vacancy clusters (VC) in zirconium. The main purpose is to produce a whole matrix of values, from small to
0168-583X/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.nimb.2004.10.052
246
D. Kulikov et al. / Nucl. Instr. and Meth. in Phys. Res. B 228 (2005) 245–249
large cluster sizes, ready for use in KMC or RT models.
2. Computational procedure The procedure for the energy calculation can be divided into three parts: (i) search for the lowest energy configuration on a rigid lattice for the cluster of selected size, defined by number of elements (vacancies) in it; (ii) relaxation and determination of the formation energy of the cluster; (iii) determination of the binding energy of a vacancy to the cluster. These three steps were automated in a shell script, which can carry out the whole calculation for clusters containing iteratively growing numbers of elements, from two to several hundreds, in a loop procedure. The calculations were performed using an empirical interatomic potential designed by Ackland and co-workers for hcp Zr [3]. 2.1. Configuration sampling We define a cluster as an object where each element has at least another element as first neighbour. Each cluster is defined primarily by its size. In order to find the most stable configuration of a selected cluster, all the different possible arrangements for the same number of elements must be considered and the corresponding total energy of the system evaluated, according to the empirical interatomic potential. The lowest energy configuration will thereby be identified and used to calculate the corresponding formation energy. A Metropolis Monte Carlo (MMC) algorithm [4] is suitable for this purpose, provided that it converges efficiently toward a total energy minimum. To speed up the calculation and enhance the probability of convergence, the possible cluster configurations at constant size were explored on a rigid lattice, i.e. by applying only atom-vacancy exchanges and computing unrelaxed total energies. It is thus postulated that a configuration exhibiting the lowest unrelaxed energy of a set of configurations of the same size will also have the lowest relaxed energy. This assumption is acceptable, considering the small distortion field associated to vacancies, provided that voids do not collapse
to loops. The limits and validity of this hypothesis are further discussed elsewhere [5]. Since full convergence cannot be always guaranteed, two different approaches were followed. The same search was repeated several times, using different initial seeds for random number generation and prolonging the simulations up to several million MMC steps. Eventually, only the search that provided the lowest configuration energy was retained. This procedure, denoted henceforth as Ôlong-searchÕ, is however too computationally demanding to be repeated systematically for a whole range of sizes and was therefore applied only to a limited number of selected cases. Alternatively, less accurate MMC searches were launched in a whole range of cluster sizes (e.g. for NV < 100): only one seed per search, all possible sizes in the considered range, all simulations stopped after the same amount of MMC steps, namely 2 · 106. In this way, although not all searches always converged, an average curve, drawn through the points corresponding to all possible clusters, can be still assumed to give a reasonable description of the formation (and binding) energies as a function of NV. This second approach will be henceforth denoted as Ôfast-searchÕ. Repeating twice or three times a fast-search with different seeds allowed in some cases to improve the results. All MMC simulations were performed in boxes containing up to 10 240 atoms (for the largest clusters), at 600 K. The CPU time required for a single long search (2–4 different seeds up to 1 · 106 MMC steps) was about 500–1500 s for one cluster, depending on its size. 2.2. Formation energy calculation The second step of the procedure is to find the formation energy of the configurations selected employing the above described methods. The most stable (lowest unrelaxed energy) configuration for each size from the previous step is therefore taken and inserted into the DYMOKA molecular dynamics (MD) code [6], where the same interatomic potential as in the previous step is used for force calculation. A block of larger size was used in this case, to minimise box size effects, which are absent on a rigid lattice with periodic
D. Kulikov et al. / Nucl. Instr. and Meth. in Phys. Res. B 228 (2005) 245–249
60
2/3
Ef = 2.72*NV
40
E f , eV
boundary conditions, but may influence the result of the relaxation. In order to relax the defective crystal, the system was quenched by forcing the velocity of each atom to zero whenever the scalar product of the velocity and the acceleration is negative. The energy of the defective system is hence calculated in the lowest energy state at 0 K and the formation energies of vacancy clusters are obtained using the following formula:
247
20
Ef ðN V Þ ¼ ðN 0 N V Þ ½Ec ðN V in ZrÞ Ec ðZrÞ; ð1Þ 0
where N0 is the total amount of atomic lattice sites in the box; Ec(NV in Zr) is the relaxed energy per atom of a Zr matrix containing a cluster of NV vacancies in its lowest energy configuration, and Ec(Zr) is the cohesive energy of pure hcp Zr and these magnitudes are the customary output of any MD code.
0
Fig. 1. The formation energies of three-dimensional VC in hcp Zr versus cluster size: raw data and interpolation.
2,5
2,0
Eb ðVÞ ¼ Ef ðclusterÞ þ Ef ðvacancyÞ
1,0
0,5
ð2Þ
where Ef (cluster) is the formation energy of a cluster of certain size, Ef (vacancy) is the formation energy of one vacancy and Ef (cluster + vacancy) is the formation energy of the cluster plus one vacancy.
3. Results and discussion The dependence of the formation energy on VC size in hcp Zr is presented in Fig. 1. The raw values can be interpolated with the following power law (checked up to a size of 1000 vacancies): Ef ðN V Þ ¼
1,5
E b, eV
Once the matrix of formation energies is known, the binding energy of a vacancy to a cluster of certain size can be obtained according to the following expression:
2=3 2:72N V :
100
NV
2.3. Binding energy calculation
Ef ðcluster þ vacancyÞ;
50
ð3Þ
One thus finds that the formation energy scales with the cluster surface.
0,0 0
10
20
30
40
50
60
70
80
90
100
110
NV
Fig. 2. Binding energies of a vacancy to a three-dimensional VC in hcp Zr versus cluster size: raw data and interpolation.
In Fig. 2 we present the values of the corresponding binding energies. The line was obtained from expression (3) using Eq. (2) and has the following expression: 2=3
EVb ðN V Þ ¼ 1:78 þ 2:72 ½N V ðN V þ 1Þ
2=3
:
ð4Þ
In this relation, the binding energy appears to be a function of the cluster radius. Fig. 2 displays a rather large dispersion around the trend curve and an analysis of the residuals (the differences between the raw data and Eq. (4) for the binding
248
D. Kulikov et al. / Nucl. Instr. and Meth. in Phys. Res. B 228 (2005) 245–249
energies was therefore made. An increasing absolute value of the residuals is observed with cluster size, with saturation for NV>30. The absolute values of the residuals as a function of cluster size (error bar) could then be approximated as 2=3
DEVb ðN V Þ ¼ 0:45 0:75 ½ðN V þ 1Þ2=3 N V :
ð5Þ
In Fig. 2, a characteristic binding energy band is noticed around 1.1 eV. The analysis of the corresponding clusters revealed that this band is related to configurations where the removed vacancy has two nearest neighbour vacancies in the same basal plane and two nearest neighbour vacancies in one of the next basal planes (see Fig. 3). The overall cluster binding energies of VC in Zr, defined as Eba ðV Þ ¼ ðEf ðclusterÞ N V Ef ðvacancyÞÞ=N V : ð6Þ This definition measures the difference between the configuration energy of the system when all vacancies form a cluster and when they are distributed in the Zr matrix without any mutual interaction. Definition (6) was used to compare our results with those reported in [7], where only small clusters were investigated and an EAM potential was used. The results of the comparison are presented in
Fig. 3. Configuration of a cluster of vacancies in which the added vacancy has a binding energy of 1.1 eV. The ‘‘bonds’’ are of course fictitious and correspond to first nearest neighbour distances between vacancies on a rigid hcp lattice.
Table 1 Comparison of overall binding energies (eV) of VC in Zr, as calculated in the present work and in [7] NV
Eba, this work (eV)
Eba, [7] (eV)
2 3 4 5 7
0.10 0.20 0.31 0.39 0.49
0.12 0.23 0.33 0.40 0.39
Table 1 and it can be seen that the values are very close. Experimentally, in irradiated Zr or Zr alloys vacancy clusters are found in the form of dislocation loop [8,9]. This is not what we find with the MMC algorithm, for limited numbers of vacancies. MMC indeed predicts the 3D void as the most stable configuration of vacancies in clusters. Further investigations are currently in course.
4. Summary and conclusions A procedure has been developed for the systematic calculation of the formation energies of vacancy clusters and corresponding binding energies in hcp Zr. This procedure is based on the use of the MMC technique on rigid lattice to identify lowest energy configurations, followed by 0 K relaxation by means of MD. This procedure neglects the role of stress in the energy minimization. However, it is acceptable for non-collapsing clusters. Systematically repeated for different cluster sizes (fast search), the proposed procedure allowed a whole matrix of formation and binding energy values to be built, up to the typical cluster diameters of interest in the case of irradiated materials. This automatic procedure caused an unphysical scatter in the corresponding binding energies. To overcome this problem, an interpolated expression for the formation and binding energies was deduced from the data. This interpolation neglects the possible existence of magic numbers for particularly compact or symmetrical cluster configurations. The interpolated formulae are of easy implementation in KMC or RE models.
D. Kulikov et al. / Nucl. Instr. and Meth. in Phys. Res. B 228 (2005) 245–249
Acknowledgements Fruitful discussions with C.S. Becquart, C. Domain and A. Legris are acknowledged. This work was performed in the framework of the convention between SCK CEN and ULB and partly under contract FIKS-CT-2001-00137 with the European communities. References [1] M.-J. Caturla, N. Soneda, E. Alonso, B.D. Wirth, T. Dı´az de la Rubia, J.M. Perlado, J. Nucl. Mater. 276 (2000) 13.
249
[2] A. Hardouin Duparc, C. Moingeon, N. Smetniansky-deGrande, A. Barbu, J. Nucl. Mater. 302 (2002) 143. [3] G.J. Ackland, S.J. Wooding, D.J. Bacon, Philos. Mag. A 71 (1995) 553. [4] O. Khrushcheva, E.E. Zhurkin, L. Malerba, C.S. Becquart, C. Domain, M. Hou, Nucl. Instr. and Meth. B 202 (2003) 68. [5] D. Kulikov, L. Malerba, M. Hou, SCK CEN Report, BLG-956, May 2004. [6] C.S. Becquart, K.M. Decker, C. Domain, J. Ruste, Y. Souffez, J.C. Turbatte, J.C. Van Duysen, Radiat. Eff. Def. Sol. 142 (1997) 9. [7] R.C. Pasianot, A.M. Monti, J. Nucl. Mater. 264 (1999) 198. [8] M. Griffiths, R.W. Gilbert, V. Fidleris, R.P. Tucker, R.B. Adamson, J. Nucl. Mater. 150 (1987) 159. [9] M. Griffith, J. Nucl. Mater. 159 (1988) 190.