Density functional theory-based derivation of an interatomic pair potential for helium impurities in nickel

Journal of Nuclear Materials 479 (2016) 240e248

Contents lists available at ScienceDirect

Journal of Nuclear Materials journal homepage: www.elsevier.com/locate/jnucmat

Density functional theory-based derivation of an interatomic pair potential for helium impurities in nickel E. Torres*, J. Pencer, D.D. Radford Canadian Nuclear Laboratories, Chalk River Laboratories, Chalk River, ON, K0J 1J0, Canada

a r t i c l e i n f o

a b s t r a c t

Article history: Received 29 January 2016 Received in revised form 23 June 2016 Accepted 6 July 2016 Available online 9 July 2016

Helium is formed in nickel as a by-product of neutron irradiation. Although helium is chemically inert and essentially insoluble in metals, under specific conditions it is known to cause metal embrittlement. Early experimental and theoretical studies on helium diffusion mechanisms have been a source of controversy. Recent density functional theory (DFT) studies of helium impurities in nickel contradict earlier theoretical studies. In this paper, a new functional form and parameters for a helium-nickel interatomic potential are proposed. The new potential used in molecular dynamics (MD) simulations correctly reproduces the relative stability of helium defects and the interstitial migration of helium in nickel. Furthermore, the computed activation energy for diffusion of helium in nickel corroborates experimental findings. The transferability of the potential is verified through a comparison with DFT predictions of the formation energies of the most stable He clusters in a Ni monovacancy. Crown Copyright © 2016 Published by Elsevier B.V. All rights reserved.

Keywords: DFT Magnetic interactions Interatomic potential Helium defects Diffusion

1. Introduction Helium is a by-product of neutron irradiation of metals in both nuclear fission and fusion reactors. In thermal spectrum fission reactors nickel is the primary source of helium, via the thermal neutron induced two step reaction 58Ni(n,g) and 59Ni(n,a). In fast spectrum fission reactors and in fusion reactors helium is produced by fast neutron irradiation through (n, a) transmutation reactions with nickel, chromium, titanium, iron and other metals. The production of helium and its aggregation into bubbles causes metal embrittlement, which has been the subject of extensive study (see Refs. [1], and references therein). Helium embrittlement of nickel alloys is of interest because some Ni-alloy components in thermal spectrum fission reactors are subjected to high neutron exposures during their long residence time within the core [2]. In addition, austenitic stainless steels and nickel-based alloys have been identified as candidates for in-core structural materials in the Generation-IV supercritical water-cooled reactor (SCWR), where they can be subject to high thermal neutron fluxes, resulting in helium production from nickel [3]. While considerable experimental effort has been devoted to understand the behavior of helium in nickel alloys and austenitic

* Corresponding author. E-mail address: [email protected] (E. Torres). http://dx.doi.org/10.1016/j.jnucmat.2016.07.009 0022-3115/Crown Copyright © 2016 Published by Elsevier B.V. All rights reserved.

steels, there have been relatively few theoretical studies in comparison to more common studies of helium in alpha-iron [1]. Early quantum mechanical electronic structure investigations were performed by Melius et al. [4] to study defects and activation energies for diffusion of helium in nickel. These studies were based on the lattice-defect method, using the Hartree-Fock approximation. The initial studies of Melius et al. formed the basis for a number of later molecular dynamics (MD) studies of helium diffusion [5] and helium bubble formation [6,7]. More recent density functional theory (DFT) [8,9] studies have been performed to determine formation energies of helium defects in nickel [10e12], and dilute austenitic Fe-Cr-Ni alloys [11]. The results of these ab initio calculations are consistent with each other, but contradict the earlier results of Melius et al. [4]. It has been found that the minimum energy interstitial location for helium is the tetrahedral position, rather than the octahedral, as previously believed, and that the barrier to interstitial helium diffusion is much lower than previously determined. It was also found that He migration between neighbor tetrahedral sites passes through an octahedral site, with a migration energy barrier of 0.129 eV [11]. The initial nucleation of helium clusters plays an important role in the formation of He bubbles. However, due the relative large size of supercell needed to appropriately model helium clusters in nickel, Empirical potential simulations are needed. As such, it is important to ensure the accuracy of the interatomic potential used for helium-

E. Torres et al. / Journal of Nuclear Materials 479 (2016) 240e248

nickel interactions in MD simulations. Despite the shortcomings of earlier simulations studies of helium in nickel, only very recently (in fact nearly coinciding with the work presented here), has there been some progress to derive a potential for helium-nickel appropriate for use in MD studies [13]. The potential of Zhang et al. is a piece-wise function derived following a method analogous to Juslin et al. [14], and was fitted in combination with the embedded-atom method (EAM) potential for nickel developed by Sheng et al. [15]. As such, while suitable for simulations of helium in pure nickel, Zhang et al. potential may be difficult to generalize to ternary mixtures of nickel, iron and chromium, which better represent high nickel alloys and austenitic stainless steels. In the present paper, we derive an accurate analytic potential for He-Ni, based on combined DFT and MD simulations. First, we investigate the relative stability and the interstitial migration barrier of helium solute impurities in pure nickel, using DFT calculations. In order to accurately model the interatomic interactions between helium and nickel atoms, we systematically analyzed the effects of spin polarization, structure relaxation, and system size on the formation energies of a helium defect in nickel. Then we report a new interatomic function, derived based on the understanding of the interactions between helium and nickel atoms, to describe He interactions in bulk nickel. The parameters of the new interatomic potential were determined by a supervised stochastic optimization algorithm using DFT results as reference data. Molecular dynamics simulations using the new He-Ni potential are shown to closely reproduce DFT formation energies of single helium defects in nickel, the energy barrier to nickel interstitial migration, and helium clusters in a single vacancy. Molecular dynamics simulations are also used with the new He-Ni potential to determine the activation energy for helium interstitial diffusion, which is found to agree with experimentally determined values.

241

Fig. 1. The volumes of the octahedrally (left panel) and tetrahedrally (right panel) coordinated nickel atoms, inside the conventional fcc unit cell, are illustrated in gray color. Interstitial sites are indicated with a solid black circle. Gray spheres represent the atoms of the nickel lattice. The pathway connecting two tetrahedral sites (red circles) is indicated by a dotted gray line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

between helium atoms and electrons were described by a nonrelativistic ultrasoft pseudopotential. We used the pseudopotentials Ni.pbe-n-rrkjus_psl.0.1.UPF and He.pbe-mt_fhi.UPF from http://www.quantum-espresso.org. In order to obtain the reference energy of a helium atom in the gas phase, a non spin-polarized calculation was performed for a helium atom including only the G-point. It was determined that a cubic supercell with a side of 6 Å, or larger, was sufficient to virtually eliminate self-interactions. In our calculations, we used a cubic supercell with a side of 26 Å to ensure sufficient accuracy. The formation energies for a single He solute impurity, at the octahedral (oct) and tetrahedral (tet) interstitial sites, (Eins), and substitutional (Esub) site (see Fig. 1), were defined as

Eins ¼ ENiþHe  Ebulk  EHe

(1)

2. Computational methods and Electronic structure calculations were performed using DFT, as implemented in the PWscf code component of the Quantum ESPRESSO integrated suite [16]. The exchange-correlation energy was described with the generalized-gradient approximation (GGA) by the Perdew-Burke-Ernzerhof (PBE) functional [17]. The valence electron wavefunction was expanded in a plane-wave basis set utilizing an upper kinetic-energy cutoff of 950 eV, which was necessary for prediction of defect formation energies within 1.0 meV of accuracy. The Brillouin zone integration was performed on an 8
8
8 Monkhorst-Pack k-point grid [18]. A first-order Methfessel-Paxton electron smearing scheme with parameter s¼0.2 eV was found appropriate for structure optimizations. Spin polarization was considered for pure nickel and with a helium solute impurity hosted. The local magnetic moments associated to different atoms were approximated by integrating the spin density within spheres centered on each atom. Total energy results were converged up to an error threshold of 0.01 meV per atom. The zeropoint energy contributions to the defect formation energies were not calculated. However, a variation of 0.01 eV has been estimated based on prior work [11]. Bulk nickel was modeled by three-dimensional periodic supercells, consisting of 2
2
2 and 3
3
3 cubic arrays of the conventional fcc unit cell (see Fig. 1), with a total of N ¼ 32 and 108 atoms, respectively. The interactions between the core of nickel atoms and valence electrons were described with a scalar relativistic ultrasoft pseudopotential, including non-linear core correction. The magnetic moments on nickel atoms were initialized to impose a ferromagnetic (fm) ordered state, but left free to relax during geometry and wave function optimizations. Interactions

Esub ¼ ENiþHe  ½ðN  1Þ=NENi  EHe ;

(2)

respectively. Here, ENiþHe is the combined total energy of bulk nickel with a helium impurity atom, ENi the energy of perfect bulk nickel, EHe the energy of an isolated helium atom, and N the total number of nickel atoms in the simulation box. The migration energy barrier (Eb ) of a helium atom moving along the pathway connecting two neighbor tetrahedral (T) sites in the [111] direction, was determined using the nudged elastic band (NEB) method. The midpoint in the [001] direction is a crowdion position, which has slightly higher formation energy than the octahedral site [11]. The middle path image was located at the octahedral (O) site; therefore the climbing image scheme was not used. The NEB calculations at the DFT level of theory were performed with the PWneb code distributed with Quantum ESPRESSO. Nudged elastic band calculations were performed in a 2
2
2 nickel supercell, with the force tolerance set to 0.1 eV/Å to control convergence. For MD simulations the Ni-Ni interactions were described using the potential of Bonny et al. [19]. This potential was chosen because it accurately reproduces defects in nickel, includes a short-range portion of the pair potential for use in cascade studies, and has cross terms enabling its use to represent nickel, iron, chromium ternary mixtures (e.g. nickel alloys and austenitic stainless steels). The He-He interactions were modeled using the Beck potential [20], which was deemed adequate for the study presented here, since it reproduces the second virial coefficient for helium, and the helium atoms in our models are either in isolation from each other

242

E. Torres et al. / Journal of Nuclear Materials 479 (2016) 240e248

or at low densities. Until the very recent potential of Zhang et al. [13], only one interatomic potential for He and Ni was described in the literature with sufficient detail to be reproduced. In the derivation of this potential, Xia et al. [7] fitted a Morse function to describe interactions using the early electronic structure calculations by Melius et al. [4]. This potential reproduces the behavior observed by Melius et al., including the octahedral position as the lowest interstitial energy for helium and ~0.6 eV as the diffusion barrier (compared to 0.43 eV found by Melius et al.). We use the potential of Xia et al. [7] for comparisons against our results obtained with the helium-nickel potential derived herein, and also include results using the recent He-Ni potential (hereafter referred to as Zhang-HeNi) by Zhang et al. [13]. The MD simulations were performed with the large-scale atomic/molecular massively parallel simulator (LAMMPS) code [21] including the GPU accelerator package [22]. Periodic boundary conditions were imposed in all directions in all the simulations. Calculations for optimization and validation were performed with supercell sizes of 3
3
3 and 6
6
6, respectively. Energies of substitutional and interstitial helium defects were calculated using the Polak-Ribìere formulation of the conjugate gradient method (as implemented in the cfg module in LAMMPS), using a force tolerance of 1.0
104 eV/Å. The energy barrier for helium interstitial migration, along the T-O-T path, was determined by the nudged elastic band method using the NEB module implemented in LAMMPS using a 6
6
6 nickel supercell. The first and last geometries, for NEB calculations, were initially optimized and stored. A total of 11 replicas were used to determine the energy barriers. The intermediate replicas were constructed by linear interpolation between first and last replica, connected by a k ¼ 10 inter-replica spring constant. The minimization was performed by a damped dynamics calculation with a force tolerance of 1.0
104 eV/Å. In order to test the accuracy of the He-Ni potential derived here, the energetics of helium clustering in a Ni vacancy were calculated with MD using nickel supercells of 20
20
20 and 2
2
2, and compared with DFT results using a 2
2
2 nickel supercell. The formation energy (Ef) of a Hen cluster in a single Ni vacancy was evaluated by the following expression

Ef ¼ ENiþHen  ½ðN  1Þ=NENi  nEHe ;

(3)

where, ENiþHe is the combined total energy of bulk nickel with a n helium atoms cluster in a single vacancy, ENi the energy of perfect bulk nickel, EHe the energy of an isolated helium atom, and N the total number of nickel atoms in the simulation box. Diffusion coefficients D for interstitial helium in nickel were also evaluated using MD, in order to make comparisons with experimental results. For diffusion simulations MD calculations were performed using a 20
20
20 supercell with a total of 32000 nickel atoms. We inserted a 1% (320) impurity concentration of helium atoms, randomly distributed in tetrahedral and octahedral sites. The system was initially relaxed by a total energy minimization, and followed by a thermal equilibration using the Langevin thermostat for 50 ps. Then the diffusion constants were obtained from the MSD computed within the canonical ensemble using the
-Hoover thermostat approximation. Up to 20 ns of total Nose simulation time was necessary for converged results. A time step 1.0 fs was used through all the simulations. In three-dimensions, the diffusion coefficient can be computed from the mean-squared displacement (MSD) of particles, using the following Einstein relation [23],

2 E D   d rðtÞ  rð0Þ 1 : D ¼ lim dt 6 t/∞

(4)

For sufficiently long simulation times, the diffusion coefficient D can be approximated as the slope of the MSD of the system as a function of time. Therefore, D can be calculated using the following equation,

MSD ¼ 6Dt:

(5)

The MSD can be directly obtained from MD simulations. The activation energy (Ea) of the migration, characteristic for diffusion to occur, can then be calculated from the dependence of the diffusion coefficient on the temperature based on the Arrhenius equation,

D ¼ D0 eEa =kB T ;

(6)

where, D0 is the pre-exponential factor, kB is the Boltzmann constant and T the temperature. Equation (6) is equivalently represented in Arrhenius plot with the following form:

lnðDÞ ¼ lnðD0 Þ  Ea ð1=kB TÞ:

(7)

Using Eq. (7), Ea and ln(D0) can be determined by linear regression of the Arrhenius plot of the diffusion coefficient obtained over a range of temperatures. 3. Results and discussion As the starting point for DFT simulations in this study, the nickel pseudopotential was tested by computing the ground state bulk properties including lattice constant (a0), cohesive energy (Ece) and magnetic moment (m). Computed bulk properties are summarized in Table 1, which includes values reported by other research groups for comparison. Here, it is seen that nickel properties obtained by us using DFT/PBE electronic structure calculations are in good agreement with DFT values reported by other research groups [10e12], based on DFT calculations performed using the planewave code VASP [25], with the Perdew and Wang (PW91) GGA exchange-correlation functional. In addition, it is noted that the DFT results are in good agreement with experimentally measured values [24]. The effects of spin-polarization, structure geometry optimization and supercell size on the formation energies of a single helium defect in nickel were also investigated. Results for various modeling parameters are summarized in Table 2, and show that the energies

Table 1 Theoretical and experimental lattice parameter (a0), cohesive energy (Ece) and magnetic moment (mB) of bulk nickel at the ground state.

DFT/PBE (this study) Experimental [24] DFT/PW91 [10] DFT/PW91 [11] DFT/PW91 [12]

a0 (Å)

Ece (eV)

m (mB)

3.52 3.52 3.519 3.522 3.52

4.84 4.44 4.84 e 4.89

0.63 0.61 0.603 0.59 0.62

Table 2 Formation energies for substitutionally and interstitially sited helium in nickel. The level of modeling, including spin-polarization, geometry relaxation and supercell size, is indicated on the top side. Energies are in eV. Size Spin Relax

2
2
2 No No

2
2
2 Yes No

2
2
2 No Yes

2
2
2 Yes Yes

3
3
3 Yes Yes

Esub Etet Eoct

4.820 7.456 7.090

3.377 6.197 5.711

4.812 6.019 6.199

3.361 4.691 4.808

3.318 4.577 4.735

E. Torres et al. / Journal of Nuclear Materials 479 (2016) 240e248

243

associated with helium in the substitutional sites are the most favorable in all cases. In order to avoid systematic errors, all the terms in the Eq. (3) were computed at the same level of theory as indicated in Table 2. From Table 2, it is seen that spin-polarization has a significant effect on the defect formation energies. The computed magnetic moment of the helium atom shows minor differences depending on the impurity site, acquiring moments of 0.0105 and 0.0024 mB on the tetrahedral and octahedral sites, respectively. A substantial quenching of the magnetic moments on the first neighbor nickel atoms, of 22% and 10%, is observed for the interstitial helium in the tetrahedral and octahedral site, respectively. In contrast, helium in the substitutional site exhibits a substantial diamagnetic behavior, with a magnetic moment 0.02 mB. However, the correct order of stability of the defects is only achieved if structure relaxation is also considered during the calculations. In fact, the percent relaxation of the volume formed by the first neighbor nickel atoms in Fig. 1, of 49.4% for the tetrahedral interstitial and 28.5% for the octahedral interstitial, indicates that local interactions play a fundamental role in the stability of helium defects. The effect of supercell size is seen by comparing the last two columns in Table 2; the relative defect formation energies do not significantly vary with the supercell size. Here, the supercells 2
2
2 and 3
3
3 have sides of 7.04 and 10.56 Å, respectively. It can therefore be concluded that possible spurious interactions between helium defect images are negligible, and that the relative stability of helium defects mainly depends on the bonding nature associated with the spatial distribution of the surrounding nickel atoms. Similar findings were reported by Zu et al., who showed that the relative stability of helium interstitial defects is primarily associated with local bonding properties between helium and neighboring nickel atoms, and does not depend significantly on the magnetic interactions [10]. The formation energies of helium defects at the substitutional and interstitial (octahedral/tetrahedral) sites computed from this study with DFT/PBE simulations are summarized in Table 3, along with corresponding DFT/PW91 results from the literature [10e12]. Table 3 also includes MD computed substitutional and interstitial energies employing the He-Ni interatomic Morse potential reported in Ref. [7] (referred hereafter only as Morse potential) and piece-wise potential from Ref. [13]. For reference, the parameters of the Morse potential are listed in Table 4. In general, there is good agreement between the DFT results. However, interstitial energies obtained with the Morse potential-based MD simulations contradict previous DFT results and those obtained in this study. The

Table 4 Parameters of the Morse function for the He-Ni interatomic pair potentials [7].rc is the pairwise cut-off radius.

Table 3 Formation energies for substitutionally (Esub) and interstitially, tetrahedral (Etet) and octahedral (Eoct) sited helium in nickel. Eb is the helium T-O-T barrier for migration. Energies are in eV.

m h i X 2 EðrÞ ¼ D0 e2aðrr0 Þ  2eaðrr0 Þ  Ai eBi r ;

Method

Esub

Etet

Eoct

Eb

DFTa DFTb DFTc DFTd MDe MDf MDg

3.318 3.205 3.185 3.23 3.736 2.722 2.908

4.577 4.393 4.46 4.5 4.501 4.532 4.499

4.735 4.556 4.589 4.65 4.616 3.928 4.653

0.116 0.17 0.129 0.15h 0.113 0.604 0.113

a b c d e f g h

DFT/PBE this work. DFT/PW91 Ref. [12]. DFT/PW91 Ref. [11]. DFT/PW91 Ref. [10]. Using the He-Ni potential from this work. Using the He-Ni Morse potential Ref. [7]. Using Zhang-HeNi potential Ref. [13]. Value calculated as Eoct  Etet.

D0 (eV)

a

r0 (Å)

rc (Å)

0.00019

1.31562

4.91572

4.5

Morse potential is able to properly reproduce the helium substitutional site as the most stable, but it is unable to reproduce the correct relative stability of interstitial defects. The recently reported piece-wise potential, while showing better agreement with DFT, underestimates the formation energy of a substitutional helium. Nevertheless, in order reduce the dependency on the available reference data, we search for a physically motivated interatomic potential function to describe He-Ni interactions, therefore, increasing the transferability of the potential to different atomic configurations. To address the significant discrepancy between molecular dynamics models and DFT, we first attempted a reparameterization of the Morse function. However, it was not feasible to obtain a satisfactory result, as it was not possible to obtain simultaneously the correct short-range asymptotic behavior of the potential and the correct order of defect energy formations. A similar conclusion was reported in Ref. [13]. Moreover, we further tested different interatomic functions, such as: Born-Mayer-Huggins, Buckingham, Beck, Gaussian, as they are implemented in the LAMMPS code [21], and their combinations. It was determined that none of these potentials, or their combinations, could be parameterized to reproduce the behavior observed in the DFT calculations. These results indicate that the interactions within the vicinity of helium atoms cannot be described by these generic functions. Thus, the effect of the contiguous nickel atoms may be underestimated. As discussed above, the relative stability of helium impurities comes primarily from effects associated with the threedimensional distribution of first neighbor nickel atoms, and not associated with the magnetic interactions, as can be observed by comparison of the 4th and 5th columns in Table 2. As such, a helium-nickel interatomic potential may be constructed without explicit representation of the magnetic interactions. However, in order to capture both the local interaction effects and asymptotic behavior, additional terms in the potential are required to adequately represent the short-range interactions. Taking this into account, we have reformulated the helium-nickel interatomic potential function. Our approach is based on a modification of the Morse potential using a linear combination of Gaussian functions, with the following analytic form:

(8)

i

where m determines the number of Gaussian functions. The fast decay of Gaussian functions makes them suitable as primitive functions to locally modify the Morse potential in the vicinity of the first nearest neighbor nickel atoms. In addition, the modified Morse potential is both continuous and differentiable everywhere, thus avoiding the problems associated with piecewise defined potentials. The LAMMPS code does not explicitly define the proposed potential function of Eq. (8), designated as Morse-mG hereafter. However, the first term, corresponding to the Morse potential, is provided. Therefore, LAMMPS was modified to include the second term of the proposed potential function, and the hybrid/overlay option available in LAMMPS was used to construct the proposed potential. In particular, m ¼ 2 was found sufficient for accurate results. The

244

E. Torres et al. / Journal of Nuclear Materials 479 (2016) 240e248

parameter space of the function expressed by the Eq. (8), for m ¼ 2, can be in general represented in vector notation as

c ¼ ðD0 ; a; r0 ; A1 ; B1 ; A2 ; B2 Þ:

(9)

The parameter optimization was performed using an iterative method developed in our group. This method is based on a random-walk algorithm in conjunction with a directed parameter search using a series of “incremental shooting” steps as a local search technique [26]. In each optimization step, the parameters were updated according to the following rule:

ciþ1 ¼ ð1 þ l*w+pÞ+ci ;

(10)

where ci is the vector of parameters at the optimization step i, l is an integer shooting parameter initialized to 1, w is a weight vector ~0.001 to 0.01, p is a random direction with components between 1 and 1, and + is the Hadamard, or vector element-wise, product [27]. The initial parameter vector c0 is provided by the user. To gauge the quality of newly generated parameters, the mean squared error (MSE) was evaluated as

 M  2 X  k k  MSEi ¼ ð1=MÞ E fci g  Eref  : 

Table 6 Total formation energy (Ef) for the most stable Hen clusters, from n ¼ 1 to 13, shown in Fig. 4. Energies are in eV. n

MDa

MDb

DFTb

1 2 3 4 5 6 7 8 9 10 11 12 13

3.73 7.21 10.53 13.77 17.27 20.60 24.25 27.75 31.35 34.85 38.19 41.52 44.70

3.75 7.30 10.77 14.25 18.03 21.77 25.94 29.96 34.24 38.47 42.53 46.89 50.98

3.36 6.38 9.40 12.35 15.72 18.88 22.58 26.33 30.08 34.94 38.93 42.59 46.39

a b

Using a 20
20
20 Ni supercell. Using a 2
2
2 Ni supercell.

80 Morse

70

Table 5 Parameters of the Morse-2G He-Ni interatomic pair potentials. (top) Morse function and (bottom) Gaussian (Gi) functions. rc is the pairwise cut-off radius. D0 (eV)

a

r0 (Å)

rc (Å)

0.00025

1.3789

4.5674

4.23

A (eV)

B

rc (Å)

15.377 7.591

0.894 0.603

4.23 4.23

G1+G2

0.10

Morse−2G

50

Zhang−HeNi

0.05

Energy (eV)

Here, M is the total number of reference data points, Ek{ci} is the MD computed energy for the ci parameters, and Ekref is the reference energy data. If MSEiþ1
G1 G2

G2

60

k

Morse

G1

(11)

40 0.00

30

2

3

4

5

20 10 0

−10 −20

0

1

2

3

4

5

Distance (Å) Fig. 2. (Solid circles) Potential energy surface of the optimized Morse-2G for He-Ni interactions. The individual Morse and Gaussian component functions are also shown. (open circles) Ni-He potential by Zhang et al. [13]. The expanded view in the inset reveals differences between the Morse-2G and Zhang-NiHe potentials in the long range region.

the Mose-2G and the Zhang-HeNi potential curves. The ZhangHeNi potential was designed for irradiation damage simulations, therefore increases very quickly for r < 1.5 Å. Our potential, on the other hand, is designed for the study of He impurities. We found that Zhang-HeNi potential, in the long range (r > 2.0 Å), does not show the correct asymptotic behavior, in fact, it becomes negative and exhibits local oscillations in this region (see inset in Fig. 2). The oscillations of the potential may result in spurious vibrations of He atoms. Nevertheless, the impact may be negligible in simulations of radiation damage involving high energy processes. In order to understand the influence of the long range interactions (see inset of Fig. 2) in the formation energies, we computed the formation energies varying the cut-off radius in the range of rc ¼ 3.0 to 6.5 Å. Fig. 5 shows the dependency of the formation energies, for the octahedral, tetrahedral and substitutional He atom, as a function of the cut-off radius. It can be clearly seen that energies computed with the Morse-2G converge for rc  4.0 Å.

E. Torres et al. / Journal of Nuclear Materials 479 (2016) 240e248

0.10

5.0 4.8

−0.10 −0.20

4.6

−0.30

4.4

−0.40

4.2

−0.50

Morse

−0.60

Fit

Energy (eV)

Energy (eV)

0.00

−0.70

Energy (eV)

0.12

Morse−2G

0.10

b

Zhang−HeNi

4.0 3.8 3.6

0.08

DFT

3.4

0.06

DFT−Fit

3.2

0.04

3.0

0.02

2.8 3.0

0.00 (T)

245

0.1

0.2

0.3

0.4

(O)

0.6

0.7

0.8

0.9

(T)

Reaction coordinates Fig. 3. NEB migration energies for a helium atom along a T-O-T path, in the [111] direction. The dotted lines are nonlinear least-square fittings to Gaussian functions.

In contrast, energies computed with the Zhang-HeNi potential substantially change with the cut-off radius up to rc  6.0 Å, and the variation in energy with the cutoff is particular prominent for the formation energy of the substitutional defect. This significant dependency with the cut-off radius is related with the oscillation of the Zhang-HeNi potential, as shown in the inset of Fig. 2. In molecular dynamics, the computational cost, time scale and the system size strongly dependent on the cut-off radius. Therefore, a small cut-off radius is a crucial characteristic of an empirical potential. The defect formation energies listed Table 3 were calculated using the Morse-2G potential. The energies in Table 3 show a good agreement between the new Morse-2G potential-based MD and

3.5

4.0

4.5

5.0

5.5

6.0

6.5

Cut−off radius (Å) Fig. 5. Formation energies of octahedral (solid diamonds), tetrahedral (solid triangles) and substitutional (solid circles) sited He atom. Energies computed with the Morse-2G and Zhang-HeNi potential [13] are indicated with solid and dotted lines, respectively.

DFT results. This result corroborates our assumption, that the Gaussian functions in conjunction with the Morse potential provide the means to obtain both, a reasonable value for the substitutional energy and the correct relative stability of interstitial energies. Nevertheless, to get the correct relative interstitial energies, a compromise in the accuracy of the substitutional energy was necessary. In Fig. 3, the molecular dynamics NEB-based migration barrier energies, computed with LAMMPS, are shown for a helium atom along a tetrahedral-octahedral-tetrahedral (T-O-T) pathway. As discussed in Hepburn et al. [11], the migration energy of interstitial helium in nickel corresponds to the difference between the formation energies of octahedral and tetrahedral interstitial helium. The migration barrier energy calculated using the Morse potential,

Fig. 4. Most stable found configurations of Hen clusters in a Ni monovacancy ranging from two (a) to thirteen (l) helium atoms. Yellow spheres represent He atoms. The first neighboring Ni atoms are represented by red spheres. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

E. Torres et al. / Journal of Nuclear Materials 479 (2016) 240e248

Fig. 3 (top), is reversed, compared to that obtained from DFT, Fig. 3 (bottom). Based on this result, it is clear that the Morse potential is not appropriate for use in atomistic simulations of helium behavior in nickel. In contrast, the migration energy computed with the Morse-2G and Zhang-HeNi potentials, Fig. 3 (bottom), shows good agreement with DFT results, with only small relative differences in the region between the interstitial sites. Furthermore, the MD energy barrier computed with the Morse-2G and Zhang-HeNi potentials, at the octahedral site, is 0.113 eV, which is in close agreement with our result of 0.116 using DFT, and is also consistent with the value of 0.129 eV found by Hepburn et al. [11], though less
table et al. [12]. Interestingly, Zhang than value obtained by Conne et al. [13] find a migration barrier of 0.12 eV using the NEB method with MD, consistent with our value, but inconsistent with the value of 0.150 eV calculated from the difference between the octahedral and tetrahedral interstitial helium formation energies, found using MD by the same group. As a further test, the lower energy Hen clusters structures in a single Ni vacancy predicted using the Morse-2G potential were compared to DFT results. A two step MD-based approach, employing the Morse-2G potential, was implemented to find the most stable Hen clusters in a single Ni vacancy. First, to construct the initial structures, a vacancy was created in the center of a 20
20
20 supercell of perfect bulk nickel, subsequently He atoms were randomly inserted within a sphere with radius of 2.4 Å centered at the Ni vacancy. In the second step, an initial MD energy minimization was performed for each cluster structure with the Ni atoms fixed, after convergence it was followed a full structure relaxation. For a comprehensive search, a set of 1000 different randomly seeded configurations were created for each cluster size. The formation energy for each set was checked for global convergence. The most stable Hen cluster structures were determined as the lowest energy structure of each set. We found that Hen cluster structures with up to n ¼ 13 maximum of helium atoms can be accommodated in a single Ni vacancy. Fig. 4 shows side view stick-ball representations of the twelve most stable configurations of Hen clusters. With the exception of He5 for which we report a new more stable geometry, some of the small clusters we found, from n ¼ 2 to 6, have similar structures to those reported by Adams et al. [6]. Nevertheless, we found slightly different geometries and orientations. On the other hand, for larger clusters (n > 6) a complete new set of structures is reported in this work. The large underestimation in the formation energies for large clusters, reported by Adams et al., can be the reason for the different findings. For comparison, the formation energies of the most stable Hen clusters were also computed at the DFT/PBE level of approximation. In these calculations a 2
2
2 nickel supercell was used. The MD and DFT total formation energies (Ef) are summarized in Table 6. The plots of the formation energy as a function of the cluster size are shown in Fig. 6. It can be observed that MD and DFT formation energies, are in close agreement, and slightly change with the supercell size. On the other hand, both approaches show an almost linear growth of the formation energy with the Hen cluster size. Nevertheless, the DFT plot shows a variation between n ¼ 9 and 10 due the magnetic effects. We observed a continuous substantial quenching of the magnetic moment of Ni atoms with the number of helium atoms from n ¼ 2 through 9. In particular, the total magnetization of nickel nearly falls to zero when going from n ¼ 9 to 10, and became zero for larger clusters. In general, we found that helium atoms inside clusters are more favorable than helium interstitials. This clearly indicates that a small helium cluster may rapidly nucleate at Ni vacancies. Molecular dynamics simulations of interstitial helium in perfect nickel showed an initial aggregation of interstitial He atoms in

50

MD (20x20x20)

45

MD (2x2x2)

40

DFT (2x2x2)

35 Ef (eV)

246

30 25 20 15 10 5 0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Helium cluster size (n)

Fig. 6. Formation energy (Ef) of Hen cluster, computed with MD and DFT, as a function of the cluster size. The supercell size used in each calculation is indicated. Energies are in eV.

small clusters. We observed the creation of interstitial Ni atoms nearby the He clusters, once clusters of n~10 of interstitial helium atoms are formed. Subsequently, He clusters take the substitutional site and further grow to larger sizes, as can be seen from Fig. 7. The created Ni interstitials are mainly located between He clusters, and are strongly correlated to the motion of neighboring helium clusters. As a final test of our developed He-Ni potential, the MD computed activation energy for thermal diffusion is compared to values reported in experimental findings. In Fig. 8 the Arrhenius plot is shown of helium diffusion constants evaluated from 400 to 700 K, in steps of 100 K, as a function of temperature. As can be seen from Fig. 8, the Arrhenius plot fit to a straight line indicates a single rate-limited thermally activated process. We determined the preexponential factor D0 ¼ 5.059
105 cm2/s and activation energy of Ea ¼ 0.20 eV from the linear regression to the Arrhenius plot. This result is remarkably close to experimental values of ~0.20 eV and

Fig. 7. He clusters (blue) of different sizes formed in nickel (red) during a MD simulation at T ¼ 400 K. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

E. Torres et al. / Journal of Nuclear Materials 479 (2016) 240e248

adjustment of the short range interactions. Finally, we have parameterized a Morse-2G helium-nickel interatomic potential to be used in conjunction with the nickelnickel potential by Bonny et al. [19], and demonstrated the applicability to the study of structures and migration of helium defects in nickel. Molecular dynamics results using the newly developed He-Ni interatomic potential are in close agreement with DFT and experimental findings. The simple analytic expression of the Morse-2G was implemented in LAMMPS using the GPU package, enabling efficient large-scale modeling of helium behavior in nickel using MD simulations. Moreover, the proposed Morse-mG function paves the way for the development of interatomic potentials of helium impurities in other metallic systems.

−13.0

−13.5

−14.0

ln(D)

247

−14.5

−15.0

−15.5

Acknowledgments −16.0 16

18

20

22

24

26

28

30

1/(kbT) (eV) Fig. 8. Plot of the logarithm of the diffusion coefficients ln(D) of helium in nickel versus 1/kBT. The reference-computed diffusion points are marked with solid diamonds.

~0.14 ± 0.03 eV activation energies for migration reported by Thomas et al. [28] and Philipps et al. [29], respectively. In contrast, earlier MD simulations reported a wide range of results, 0.48 [7], 0.53 [30] and 0.66 eV [5], despite being based on interatomic potentials all derived from the early results by Melius et al. [4]. 4. Summary and conclusions In this paper, a systematic DFT study of the effects of spin polarization, structure relaxation, and system size on the formation energies of single helium defects in nickel was presented. It was shown that helium defect formation energies in nickel are mainly determined by short-range interactions between helium and the first neighboring nickel atoms. In particular, we found that although magnetic interactions have an important effect in lowering the formation energies, they are not the determining factor of the relative stability of helium impurities in nickel. A previously reported He-Ni Morse potential, used in molecular dynamics studies such as Ref. [7], was examined and it was shown that it cannot reproduce the relative stability of helium defects, and consequently, the barrier for interstitial helium migration is in reverse. Moreover, the Morse-based prediction of stability for a substitutional helium defect (i.e. energy barriers between substitutional and interstitial positions) is considerably higher than that predicted by the ab initio results. A more recently devised He-Ni pair potential parameterized by Zhang et al. [13] presents improvements, however, shows inconsistencies in the interstitial migration energy. Furthermore, oscillations below zero of this potential may results in spurious vibrations of He atoms in nickel. Notwithstanding, for high energy studies this may not be an issue. In order to better represent the local interactions contributing to the relative stability of interstitials and diffusion barrier, we attempted a reparameterization of the Morse function. However, it was not possible to reproduce simultaneously the correct relative stability of the helium defects and the correct asymptotic behavior. Since the relative formation energies of interstitial helium impurities are mainly dominated by local interactions between helium and the first nearest neighbor nickel atoms, we have proposed a new potential: The Morse-mG potential, expressed by a Morse function superimposed with a linear combination of Gaussian functions. This new function is able to reproduce both short and longer range interactions via the Gaussian terms, which enable

We would like to thank Dr. Trung Nguyen for his assistance with the GPU accelerated implementation of the new potential function in LAMMPS. The authors are grateful to the advanced computing facilities at the Canadian Nuclear Laboratories (CNL) for providing generous support to this work. References [1] Y. Dai, G.R. Odette, T. Yamamoto, 1.06-The effects of helium in irradiated structural alloys, in: R.J.M. Konings (Ed.), Comprehensive Nuclear Materials, Elsevier, Oxford, 2012, pp. 141e193. [2] R. Boothby, 4.04-Radiation in nickel-based alloys, in: R.J.M. Konings (Ed.), Comprehensive Nuclear Materials, Elsevier, Oxford, 2012, pp. 123e150. [3] D. Guzonas, R. Novotny, Supercritical water-cooled reactor materials Summary of research and open issues, Prog. Nucl. Energy 77 (2014) 361e372. [4] C.F. Melius, C.L. Bisson, W.D. Wilson, Quantum-chemical and lattice-defect hybrid approach to the calculation of defects in metals, Phys. Rev. B 18 (4) (1978) 1647e1657. [5] A. Adams, W. Wolfer, On the diffusion mechanisms of helium in nickel, J. Nucl. Mater 158 (1988) 25e29. [6] J. Adams, W. Wolfer, Formation energies of helium-void complexes in nickel, J. Nucl. Mater 166 (3) (1989) 235e242. [7] J. Xia, W. Hu, J. Yang, B. Ao, X. Wang, A comparative study of helium atom diffusion via an interstitial mechanism in nickel and palladium, Phys. Status Solidi B 243 (3) (2006) 579e583. [8] P. Hohenberg, W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136 (3B) (1964) B864eB871. [9] W. Kohn, L.J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140 (4A) (1965) A1133eA1138. [10] X.T. Zu, L. Yang, F. Gao, S.M. Peng, H.L. Heinisch, X.G. Long, R.J. Kurtz, Properties of helium defects in bcc and fcc metals investigated with density functional theory, Phys. Rev. B 80 (5) (2009) 054104. [11] D.J. Hepburn, D. Ferguson, S. Gardner, G.J. Ackland, First-principles study of helium, carbon, and nitrogen in austenite, dilute austenitic iron alloys, and nickel, Phys. Rev. B 88 (2) (2013) 024115.
table, D. Monceau Andrieu, First-principles nickel database: ener[12] D. Conne getics of impurities and defects, Comput. Mater. Sci. 101 (2015) 77e87. [13] W. Zhang, C. Wang, H. Gong, Z. Xu, C. Ren, J. Dai, P. Huai, Z. Zhu, H. Deng, W. Hu, Development of a pair potential for NiHe, J. Nucl. Mater 472 (2016) 105e109. [14] N. Juslin, K. Nordlund, Pair potential for FeHe, J. Nucl. Mater 382 (23) (2008) 143e146. [15] H. Sheng, M. Kramer, A. Cadien, T. Fujita, M. Chen, Highly optimized embedded-atom-method potentials for fourteen fcc metals, Phys. Rev. B 83 (2011) 134118. [16] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A.P. Seitsonen, A. Smogunov, P. Umari, R.M. Wentzcovitch, Quantum espresso: a modular and open-source software project for quantum simulations of materials, J. Phys. Condens. Matter 21 (39) (2009) 395502 (19pp). [17] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (18) (1996) 3865e3868. [18] H.J. Monkhorst, J.D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13 (12) (1976) 5188e5192. [19] G. Bonny, N. Castin, D. Terentyev, Interatomic potential for studying ageing under irradiation in stainless steels: the FeNiCr model alloy, Model. Simul. Mater. Sci. Eng. 21 (8) (2013) 085004. [20] D. Beck, A new interatomic potential function for helium,, Mol. Phys. 14 (4)

248

E. Torres et al. / Journal of Nuclear Materials 479 (2016) 240e248

(1968) 311e315. [21] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comp. Phys. 117 (1) (1995) 1e19. [22] W.M. Brown, S.J.P.P. Wang, A.N. Tharrington, Implementing molecular dynamics on hybrid high performance computers - short range forces, Comp. Phys. Comm. 182 (2011) 898e911. €rme [23] A. Einstein, Über die von der molekularkinetischen Theorie der Wa geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen, Ann. Phys. 322 (8) (1905) 549e560. [24] C. Kittel, Introduction to Solid State Physics, Wiley, Hoboken, NJ, 2005. [25] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (16) (1996) 11169e11186.

[26] E. Torres, G.A. DiLabio, Density-functional theory with dispersion-correcting potentials for methane: bridging the efficiency and accuracy gap between high-level wave function and classical molecular mechanics methods, J. Chem. Theory Comput. 9 (8) (2013) 3342e3349. [27] G.P.H. Styan, Hadamard products and multivariate statistical analysis, Linear Algebra Appl. 6 (1973) 217e240. [28] G.J. Thomas, W.A. Swansiger, M.I. Baskes, Low temperature helium release in nickel, J. App. Phys. 50 (11) (1979) 6942e6947. [29] V. Philipps, K. Sonnenberg, Interstitial diffusion of He in nickel, J. Nucl. Mater 114 (1) (1983) 95e97. [30] M.I. Baskes, C.F. Melius, Pair potentials for fcc metals, Phys. Rev. B 20 (8) (1979) 3197e3204.