ARTICLE pubs.acs.org/JPCC
Bond-Order Potential for Erbium-Hydride System S. M. Peng,† L. Yang,*,‡ X. G. Long,† H. H. Shen,†,‡ Q. Q. Sun,‡ X. T. Zu,‡ and F. Gao*,§ †
Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics, Mianyang 621900, China Department of Applied Physics, University of Electronic Science and Technology of China, Chengdu 610054, China § Pacific Northwest National Laboratory, P.O. Box 999, Richland, Washington 99352, United States ‡
ABSTRACT: Interatomic potentials for an ErH system are derived based on an analytical bond-order scheme. The model potentials provide a good description of the bulk properties and defect properties of hcp-Er, including lattice parameters, cohesive energy, elastic constants, point defect formation energies, and surface and stacking fault energies. In addition to experimental data, an ab initio method is used to construct the necessary database of different phases. We demonstrate that such potentials can reproduce the hydrogen behavior in an α-phase ErH system for a low hydrogen/metal ratio. Especially, the present potentials can be employed for modeling the energetics and structural properties of fcc ErH2, including lattice parameters, elastic constants, bulk modulus, Young’s modulus, shear modulus, as well as the formation energies and migration barriers of point defects in ErH2.
I. INTRODUCTION Interest in rare-earth hydrides has been stimulated in the past few years for several reasons.1,2 One of the most interesting reasons is that rare-earth metals present a large affinity for hydrogen.36 When the rare-earth metals react with H, a set of very similar phase diagrams is formed with various H concentrations ranging from 0 to 3 in atomic ratio.6 Among the rare-earth metals, Er metal has been shown to be an excellent material for accommodating H even at low temperatures,2 and the phase diagram of ErH is typical for the lanthanide family of metalH systems.7 The solidsolution (α) phase of the ErH system is generally formed at low H concentrations in the metal lattice, retaining the hexagonal structure of the host metal, but the coexistence of α-phase with β-phase (fcc structure) ErH2 is observed at higher H concentrations, followed by a narrow region of pure ErH2. A region of mixed β-phase ErH2 and hexagonal trihydride phases persists with increasing H concentration up to 2.7 in atomic ratio, above which a pure hexagonal phase of ErH3 is observed.8 Recently, first principles methods were used to calculate the site occupation and migration of H in β-phase ErH29 and in α-phase ErH.10 The dynamics simulations of H in Er at the atomic scale are useful for understanding the behavior of H in an ErH system. It is well-known that the reliability of the simulations greatly depends on the accuracy of the interatomic potential used in the simulations. Modified embedded atom method (MEAM) potentials for ErEr interaction have been developed by Baskes et al.11 and Hu,12 but there is no interatomic potential describing the ErH interaction in the literature. The ErEr MEAM potential proposed by Baskes et al. was fitted to several physical parameters, including lattice constants, elastic constants, cohesive energy, vacancy formation energy and the bcc-hcp structural energy difference, reproducing this extensive database, r 2011 American Chemical Society
but the properties of point defects were not well reproduced. For example, the binding energy of a divacancy was found to be negligible in many hcp metals, including Er, whereas the ab initio calculations in the present work show that a divacancy in Er is bound with a positive binding energy. For Hu’s potential, the thermodynamic properties of MgEr alloys were considered, but the interatomic potential of pure Er was not provided. Therefore, it is difficult to evaluate Hu’s potential and its application to study defect properties in Er. Juslin et al.13 have developed a HH interaction based on an analytical bond-order potential (BOP). To fully take advantage of their HH potential, we have constructed the ErEr and ErH interatomic potential based on the same analytical BOP scheme in the present work. The BOP approach is based on the TersoffBrenner potential1416 that has been used to describe not only SiC, WCH, WH, BeCH, PtC, and ZnO1721 but also transition metals, such as Mo, Fe, and W.2224
II. METHOD A. Potential Formalism. The general formula of the analytical bond-order potential is given in detail in refs 1621. In this approach, the total energy of a system is written as a sum over individual bond energies by bij þ bji A c R Vij ðrij Þ fij ðrij Þ Vij ð1Þ E¼ 2 i>j
∑
Received: September 19, 2011 Revised: November 3, 2011 Published: November 04, 2011 25097
dx.doi.org/10.1021/jp2090523 | J. Phys. Chem. C 2011, 115, 25097–25104
The Journal of Physical Chemistry C
ARTICLE
Table 1. Computational Cell Size and Number of k-Points in the Present PW91/DFT Calculations for Different Properties, where “Gamma” and “MonkhorstPack” Represent the k-Mesh Generation Methods in the Brillouin Zone property
computational cell size 15 Å 15 Å 15 Å
k-points
dimer
ErEr
elastic constant
ErH Er
1 (Gamma)
One cell (2 atoms)
11 11 9 (Gamma)
ErH2
One cell (12 atoms)
9 9 9 (Monkhorst-Pack)
defect properties in Hcp-Er and ErH
4 4 3 cells (96 atoms)
2 2 2 (Gamma)
surface energy of Er
9 layers
9 9 1 (Gamma)
2 2 atoms per layer 15 Å vacuum stacking fault energy of Er
12 layers (ABABABCBCBCB)
I1
9 9 1 (Gamma)
2 2 atoms per layer 15 Å vacuum 12 layers (ABABABCACACA)
I2
9 9 1 (Gamma)
2 2 atoms per layer 15 Å vacuum 2 2 2 cells (96 atoms)
properties in Fcc-ErH2
The pair-like repulsive and attractive terms are taken as Morselike pair potentials pffiffiffiffiffi D0 exp½ β 2Sðr r0 Þ V R ðrÞ ¼ ð2Þ S1 VA ¼
pffiffiffiffiffiffiffi SD0 exp½ β 2=Sðr r0 Þ S1
ð3Þ
where D0 is the dimer bond energy, r0 is the dimer bond distance, and S and β are adjustable parameters. The interaction range is restricted by a cutoff function 8 > 1, > r e RD > < 1 1 π sin ðr RÞ=D , jR rj e D ð4Þ f c ðrÞ ¼ 2 2 2 > > > : 0, rgR þ D where D and R are adjustable parameters and determine the cutoff range and interval. The bond-order parameters, including three-body contributions and angularity, are described by bij ¼ ð1 þ χij Þ1=2 χij ¼
fikc ðrik Þgik ðθijk Þ exp½αijk ðrij rik Þ ∑ kð6¼ i, jÞ
ð5Þ ð6Þ
where αijk is an adjustable parameter and the sum is over the third atom k. The angular function g(θ) has the form " # c2 c2 gðθÞ ¼ γ 1 þ 2 2 ð7Þ d d þ ðh þ cos θÞ where γ, c, d, and h are adjustable parameters. B. Constructing Fitting Database. The fitting database is generally constructed from the available experimental data for hcp-Er. However, experimental data for ErH systems are very scarce due to its chemical instability. As described above, the hexagonal structure ErH system is formed at low H concentrations, but the coexistence of hexagonal ErH with fluorite structure ErH2 is observed at higher H concentrations, followed by a narrow region of pure ErH2. To ensure that a potential is transferable for different
4 4 4 (Monkhorst-Pack)
phases, the fitting database has to cover the possible structures of different atomic coordinations and different phases. For these reasons, ab initio calculations based on the density functional theory (DFT) framework have been carried out to determine several Er and ErH structures. The database was constructed to include cohesive energies, lattice constants, elastic constants, bulk moduli and defect formation energies for Er and ErH2, as well as Young’s moduli and shear moduli of ErH2. The PW91/DFT calculations were carried out using the Vienna ab initio simulation package (VASP) code.25,26 The projector-augmented wave (PAW) method was used to determine the interactions between ions and electrons,27,28 while exchange and correlation were described by the generalized gradient approximation (GGA) of Perdew and Wang (PW91).29 The standard pseudopotentials of Er and H were taken from the VASP library, which provides a good description of the properties of a perfect crystal. The plane-wave basis sets with energy cutoff of 500 eV were performed throughout the present work. During relaxations, the Brillouin zone integration was achieved using a MethfesselPazton smearing of σ = 0.1 eV. In the ab initio calculations, the results significantly depend on the computational cell size and the number of k-points that are displayed in Table 1 for the present PW91/DFT calculations in detail. The calculations for an Er dimer and an ErH dimer were performed using a supercell with a border length of 15 Å. In the calculations the two atoms of a dimer were first separated at a distance in the supercell, and then relaxed using the conjugate gradient algorithm until the total energy of the dimer and the force on every atom were converged to 1 105 and 1 103 eV, respectively. In addition, spin polarized calculations were performed. In the calculations of elastic constants, Gamma-centered grids with 11 11 9 k-points were performed for Er, whereas the MonkhorstPack scheme30 with 9 9 9 k-points was used for ErH2. In the calculations of defect properties in hcp Er and ErH systems, Gamma-centered grids with 2 2 2 k-point meshes were used in a 96 atom supercell. However, the 4 4 4 k-point meshes were used to check the coverage of defect calculations, and it is found that the difference in defect formation energy is less than 0.02 eV. To calculate the surface energies of Er crystals, a slab with nine (2 2) Er atomic layers, with a 15 Å vacuum between any two 25098
dx.doi.org/10.1021/jp2090523 |J. Phys. Chem. C 2011, 115, 25097–25104
The Journal of Physical Chemistry C
ARTICLE
successive metal slabs, was imposed. Gamma-centered grids with 9 9 1 k-points within the surface Brillouin zone was used. The surface energy is determined by Efsurf ¼
ETsurf NEc 2A
parameter
ð8Þ
where ETsurf is the total energy of a computational cell with two free surfaces, N is the number of atoms, Ec is the cohesive energy, and A is the surface area. During the calculation of stacking energy in Er crystals, a slab with twelve (2 2) Er atomic layers and a space of 15 Å between periodically repeated slabs were used to void the self-interaction between the stacking faults.31 The stacking sequence ABABABCBCBCB is employed for I1 and ABABABCACACA for I2. The stacking energy is defined as Efstack ¼
ETstack NEc 2Efsurf ðbasalÞ A
ð9Þ
where ETstack is the total energy of a computational cell with a stacking fault, and Efsurf(basal) is the surface energy of the basal plane. N and Ec are the same as those in eq 8. A is the stacking fault area. To determine the structure of ErH2, Brillouin zone (BZ) sampling was performed using the MonkhorstPack scheme30 with 4 4 4 k-point meshes in a 96 atom supercell. In addition, the climbing image nudged elastic band (CI-NEB) method32 was used to investigate the migration mechanisms of defects and H in a supercell with fixed volume and shape. C. Fitting Methodology. Once an extensive fitting database was constructed, the methodology for determining the set of potential parameters that reproduce as many properties as possible was initiated. Starting from an initial guess, the potential parameters are optimized by adjusting the properties predicted by the potential to fit the experimental data, including the lattice parameters, elastic constants, cohesive energy and defect formation energies, using a conjugate gradient least-squares minimization algorithm. The objective function, U, which determines if each individual set of fitting parameters can be accepted, is defined as U ¼
∑i wi ½fi ðλnÞ Fi 2
Table 2. Potential Parameters for ErEr and ErH Interactions in the Present Work, along with Those for the HH Interaction Given by Brenner15a
ð10Þ
where fi represents the calculated property that depends on the potential parameters λn. The weight of each data item i in the fitting process, denoted as wi, determines how well the fitting can reproduce the corresponding property. The fitting code first produces a set of parameters λn, which can be used to calculate the fitting properties fi (λn). The MD code LAMMPS33,34 was employed to obtain the formation energies of defects and struc√ ture properties. In the fitting procedure, a 5a 4 3a 4c (where a and c are the lattice parameters of the Er crystal) cubic box containing 320 Er atoms was used for the hcp ErEr and ErH systems, and a 3a 3a 3a (where a is the lattice parameter of fcc ErH2) simulation box containing 324 atoms was used for the fcc ErH2 system. The objective function was determined until all the fitting properties were calculated. This process was repeated until the selection criteria was satisfied. When a satisfactory set was found, the potential was tested against the properties that were not included in the fitting, and if not performing adequately, as compared to experimental data or PW91/DFT calculations, a new parameter set was refitted.
ErEr
ErH
HH (Brenner)
D0 (eV)
1.227637
2.853565
4.7509
R0 (Å)
3.185783
1.834452
0.74144
β (Å1)
0.946477
1.169382
1.9436
S
2.857930
2.723064
γ c
0.151171 0.073741
0.170183 8.023344
12.33 0.0
0.127005
1.809234
1.0
d h
0.47288
2.3432
0.512327
1.0
R (Å)
5.601367
3.920701
1.40
D (Å)
0.171007
0.100477
0.30
Here αErErEr = 1.488350, αErErH = 1.218961, αEHrEr = 1.141117, αErHH = 0.169670, αHErEr = 0.550607, αHErH = 0.0, αHHEr = 0.0, and αHHH = 4.0. a
The parameters derived for the ErEr and ErH interactions are given in Table 2, where the parameters for the HH interaction from Brenner15 are also provided.
III. RESULTS A. ErEr. In order to evaluate the reliability of the above potentials, we first calculate the various physical properties of bulk hcp-Er, including lattice constants, elastic constants, structural energy differences, point defect properties (vacancy and interstitial formation energies, divacancy binding energy, vacancy migration energy), planar defect properties (surface energy and stacking fault energy), and melting point. Except for the melting point, all the calculations√in this section and the next section are carried out with a 10a 6 3a 6c cubic box containing 1440 Er atoms. These properties calculated from the present BOP are shown in Table 3, together with data from experiments and PW91/DFT calculations. Here, the properties marked with a “*” are those used in the fitting process. The corresponding data from the modified embedded atom method (MEAM)11 are also presented for comparison. As can be seen from Table 3, the cohesive energy of hcp-Er calculated from the present potential is 3.292 eV/atom, which is in excellent agreement with the experimental value (3.29 eV/ atom),35 but significantly higher than the PW91/DFT data (4.28 eV/atom). Similar to the MEAM calculations,11 the present potential well reproduces the bulk modulus and the elastic constants, and the differences of the elastic constants between the present calculations and the experimental results37 are less than 3.5%. The structural energy differences (ΔEhcpfbcc, ΔEhcpffcc, ΔEhcpfsc, and ΔEhcpfdia) obtained in the present work are also listed in Table 3, along with the values obtained from other calculations. Among other target properties, the structural energy differences between hcp and other phases were assigned relatively low weights during the fitting procedure. It is found that the fitted hcp lattice is indeed the most stable structure, although, the structural energy differences are small for fcc and bcc compared to PW91/DFT data and those of the MEAM calculations. Agreement with PW91/DFT and MEAM11 for simple cubic (sc) and diamond (dia) is excellent. In terms of point defect properties, only the vacancy formation energy is experimentally available in hcp-Er, and the present BOP 25099
dx.doi.org/10.1021/jp2090523 |J. Phys. Chem. C 2011, 115, 25097–25104
The Journal of Physical Chemistry C
ARTICLE
Table 3. Calculated Physical Properties of Er Using the Present BOP, in Comparison with Those Obtained from Experiments and Calculated by PW91/DFT Methods and Other Potentialsa present BOP
experimental
PW91/ DFT MEAM11
*Ec
3.292
3.29b
4.28
3.04
*a
3.559
3.559c
3.562
3.549
*c
5.592
5.592c
5.550
*C11
84.99
86.78d
85.66
86.56
*C12
27.79
28.51 d
25.08
27.28
*C13 *C33
24.63 82.03
24.18 d 82.05 d
23.78 90.99
24.08 85.43
*C44
25.61
26.48 d
29.98
27.23
*B
45.12
45.48
45.5
45.49
*ΔEhcpfbcc
0.04
0.15
0.207
*ΔEhcpffcc
0.002
1.95
0.054
*ΔEhcpfsc
0.50
0.83
0.46
*ΔEhcpfdia
1.27
*Efvac Em vac
1.36 0.76
1.77 0.62
1.37
in plane
5.592
1.54 1.32e
out of plane 0.77
0.65
0.18
0.23
0.68
out of plane 0.17
0.24
0.64
T
3.37
3.06
O
3.58
2.66
C
3.32
2.93
S BO
3.53 3.04
3.06 2.31
BS
3.02
2.80
BC
unstable
unstable
BT
3.49
3.37
basal
0.623
1.047
2.328
prism
0.781
1.197
2.238
Efstack
I1
0.016
0.031
dimer
I2 *Ec
0.034 1.23
0.081 1.88
*r0
0.319
Ebdivac EfI
Efsurf
Tm
in plane
1350 ( 50
Figure 1. Eight self-interstitial configurations in an hcp structure. O stands for octahedral, C for crowdion, S for split dumbbell, and BT, BO, BC, and BS for the corresponding positions on the basal plane, where the black spheres show the hcp lattice sites and the green spheres represent the vacant lattice sites in the S and BS configurations.
paths of a vacancy have been considered in the present study, the in-basal plane (in plane), and the out-of-basal plane (out of plane). The migration energy for a vacancy is the difference between the energy for an atom at the saddle point and that at its equilibrium site as it moves from its crystal site to the nearest vacant site.12 It is clear that the migration energy in plane is slightly lower than that out of plane. Similarly, two configurations of a divacancy (in plane and out of plane) have also been considered to determine their formation energies and the corresponding binding energies. The binding energy of a divacancy can be evaluated by19 Eb2vac ¼ 2Efvac Efdivac
ð11Þ
where Efvac and Efdivac are the formation energies of a single vacancy
0.160
0.311 1802f
a
The properties include the cohesive energy Ec (eV/atom), the lattice parameters of a and c (Å), the elastic constants Cij and bulk modulus B (GPa), the structural energy differences ΔE (eV/atom), the relaxed vacancy formation energy Efvac, vacancy migration energy Em vac, divacancy binding energy Ebdivac (eV), the self-interstitial formation energy EfI (eV), the basal and prism plane surface energies Efsurf (J/m2), the stacking fault energy Efstack (J/m2), the ErEr dimer bond distance r0 (nm), and the melting point Tm (K) of an hcp-Er. Properties marked with a “*” are those used for fitting in the present work. b Reference 35. c Reference 36. d Reference 37. e Reference 38. f Reference 39.
well reproduces the experimental value of 1.32 eV.38 However, the vacancy formation energy obtained with the present potential is 1.36 eV, which is 0.39 eV lower than the PW91/DFT data of 1.75 eV. This suggests that the PW91/DFT calculation overestimates the vacancy formation energy by about 0.39 eV. In order to have a self-consistent data set, the formation energies associated with vacancies have been modified, which will be discussed in section B1 of the results. In Table 3, two migration
and a divacancy, respectively. It is clear that both the migration energy of a vacancy and the binding energy of a divacancy in the out-of basal plane are similar to those in the basal plane, which agrees with the results from the PW91/DFT calculations. It is of interest to note that the present potential gives a positive binding energy for a divacancy, which indicates that the divacancy configuration is stable. However, a negative value of the binding energy was predicted in the potential model proposed by Baskes and Johnson.11 In Table 3, the formation energies of several self-interstititial atoms (SIAs) (EfI) are also evaluated, and the corresponding initial configurations of SIAs are shown in Figure 1, where T denotes a tetrahedral site, O represents an octahedral site, the C site is the midway between two nearest-neighbor atoms out of the basal plane, and S indicates a Æ0001æ split dumbbell. Basal tetrahedral (BT) and basal octahedral (BO) are the projections of T and O, respectively, onto a nearby basal plane. BS involves two atoms equidistant from a vacant lattice site, while BC is the midpoint between two nearest neighbors in the basal plane. Except for the BC site, they are all stable using the potentials developed in this work, which agrees with the PW91/DFT results. It should be noted that the most stable configuration for the SIAs appears to be BO, followed by the O and BS interstitials for the PW91/DFT calculations.10 It is clear that the formation energies of self-interstitials are larger than those obtained by the PW91/ DFT, but the BOP potential gives almost the same ground state interstitial configuration when compared to the PW91/DFT calculations, because the formation energy of BO is only 0.02 eV larger than that of the most stable structure of BS. The present potential can thus be considered to be reasonable in describing the properties of the bulk Er and point defects in hcp-Er. 25100
dx.doi.org/10.1021/jp2090523 |J. Phys. Chem. C 2011, 115, 25097–25104
The Journal of Physical Chemistry C Planar defect properties are also quite important to test an interatomic potential. We thus calculate the unrelaxed surface energy and stacking fault energy of hcp-Er using eqs 8 and 9. The surface energies of the basal and prism planes determined from the present potential are 0.623 and 0.781 J/m2, respectively. However, the surface energies for the two planes determined using the MEAM are roughly the same, with the values of 2.328 and 2.238 J/m2. As compared with the PW91/DFT values of 1.043 and 1.181 J/m2, the present potentials are more reasonable in describing the planar defect properties. The stacking faults that are related to the dislocation glide and loop-habit planes provide a further pathway to test interatomic potentials. On the basal plane, there are three possible stacking faults. There are two intrinsic faults having the stacking sequences ABABCBCB (I1) and ABABCACA (I2), and an extrinsic fault with the stacking sequence ABABCABAB (E). I2 is the most important fault in connection with plastic deformation (dislocations splitting into Shockley partials), and therefore, the most reported values of stacking fault energy are associated with this fault. In this work, the stacking fault energies of two intrinsic faults have been considered and are listed in Table 3. In the BOP fitting, the cohensive energy of dimers, i.e., the dimer bond energy D0, and the bond distance r0 of dimers are regarded as adjustable parameters, such as for the ZnZn iteraction.21 In the present work, we also considered the Er dimer properties, but among the other target properties, the D0 and r0 of an Er dimer were given relatively low weights during the fitting procedure. It is of interest to note that the final values of D0 = 1.23 eV and r0 = 0.319 nm are comparable to the PW91/DFT values of 1.88 eV and 0.311 nm for the dimer. The DFT calculations on the dimer can be found in B: Constructing fitting database of the METHOD. To understand thermal dynamic properties, the melting behavior of Er has been investigated using the present potential, which was determined by simulating a solidliquid interface at zero pressure and different temperatures.18,40 The simulations were performed with a 16 4 4 simulation box containing 1024 Er atoms, where temperature and pressure were controlled using the Berendsen thermostat and barostat.41 The left part of the constructed solidliquid interface is fixed, while the right part is relaxed at 3000 K for 5 ps, followed by cooling down and quenching the system to 0 K. The initially formed state is then relaxed at different temperatures, and equilibrated for times up to 0.5 ns. The present potential predicts a melting temperature of 1350 ( 50 K, 452 K lower than experimental value (1802 K).39 However, it should be noted that the melting temperature may be related to the formation energies of SIAs and vacancies, as discussed by Bacon42 and Hu.43 To understand this, we have refitted the potential with higher defect formation energies of Er, but the results show that the melting temperature remains almost unchanged with increasing formation energies of point defects, which is similar to the melting behavior of hcp Ru and Re.43 B. ErH. B.1. H in hcp-Er. The phase diagrams of some rareearth-H systems, RH(D)x (where R = Sc, Y, Ho, Er, Tm, and Lu), are characterized by a solid solution phase extending to rather high H concentrations, which is stable down to 0 K.44 The HEr phase diagram clearly demonstrates that when the ratio of H atoms to Er atoms (H/Er) in a HEr system drops below 1.85, the precipitation of the hexagonal α-phase is thermodynamically favorable.9 In order to understand the H behavior and its thermodynamics in α-phase metal-H systems, the ErH potential should reproduce the H properties in hcp-Er. The potential parameters are
ARTICLE
Table 4. Formation Energies (eV) of a Single H Atom in Tetrahedral (T), Octahedral (O), Basal Tetrahedral (BT), and Basal Octahedral (BO) Interstitial Sites and Substitutional (S) Site, as well as a Pair of H Atoms in an hcp-Er, in Comparison with Those from PW91/DFT Calculationsa configuration
present BOP
PW91/DFT
*T
0.99
1.00
*O *BT
0.82 0.90
0.77 0.90
*BO
0.20
0.18
S
1.84
2.05
T1-T3
2.05
2.10
T1-T4
2.01
1.96
T2-T3
1.99
2.02
T1-O
1.67
1.71
T1-T2
1.19
1.53
a
The formation energies of H point defects refer to half of the total energy of an H2.
presented in Table 1, and the properties of H defects are given in Table 4, along with the PW91/DFT data for comparison. Similar to H in hcp-Be,46 all of the substitutional (S) sites are equivalent in an hcp structure, whereas four well-known interstitial sites are potential candidates for hydrogen atom locations in hcp metals, including T, O, BO, and BT, as shown in Figure 1. The formation energies of H at the four sites as determined by the present potential are 0.99, 0.82, 0.90, and 0.20 eV, respectively, which are in good agreement with those obtained from the PW91/DFT calculations. All of the calculations refer to half of the total energy of H2. The formation energy of an interstitial H is defined to be47 EfH ¼ Etot ½H Eper ½Er εH
ð12Þ
where Etot[H] is the total energy of a computational cell containing an interstitial H, Eper[Er] is the total energy of a perfect hcp Er crystal, and εH is the energy per H atom in an H molecular. The potential gives EfT,H < EfBT,H < EfO,H < EfBO,H, which suggests that the ground state is a tetrahedral (T) position. In order to correctly determine H-vacancy formation energies, we set the potential parameters R and D large enough to ensure the correct description for the ErH interactions, particularly for the case of a substitutional H atom. From Table 4, it is clear that the formation energy of a substitutional H is in excellent agreement with the PW91/DFT data. Note that the value of the substitutional H formation energy from the PW91/DFT calculations in Table 4 is 0.39 eV lower than the original value of 2.44 eV reported in ref 10. This is due to the fact that the formation energy of a vacancy obtained with the present potential is 1.36 eV, which is 0.39 eV lower than that from the PW91/DFT calculations (1.75 eV). Because we aimed to accurately fit the binding energy of an H atom to a vacancy, rather than its substitutional formation energy, the PW91/DFT value was adjusted to consider this difference.48 Though the formation energy (1.84 eV) of a substitutional H calculated by the present potential is lower than the PW91/DFT value of 2.05 eV, both the PW91/DFT and empirical potential calculations show that the formation energy of a substitutional H atom is larger than those of the interstitial H atoms, suggesting that an interstitial solid solution of H is favorable at relatively low H concentrations in hcp-Er, which is in agreement with experimental observations.49 25101
dx.doi.org/10.1021/jp2090523 |J. Phys. Chem. C 2011, 115, 25097–25104
The Journal of Physical Chemistry C
ARTICLE
Table 5. Properties of the Bulk Phases of ErH2 from the Present BOP, in Comparison with Those from Experiments and PW91/DFT Calculationsa PW91/DFT Sandia
Sandia
this
experiment
Snow45
Wixom11
work
5.126b
5.1295
present BOP
Figure 2. Several interstitial sites in an hcp structure, where the squares and circle present T and O sites, respectively. The subscript numbers in T show different T interstitial sites.
Hydrogen ordering in many hcp metals was also observed by neutron-diffraction and diffuse-scattering, which revealed that the H atoms in pairs occupy the second-neighbor T interstitial sites along the Æ0001æ direction around a metal atom and form quasi-linear chains along this direction.50,51 In the present BOP, we also consider the properties of an H pair in bulk Er. Several possible interstitial sites are shown in Figure 2,47 i.e., two H atoms occupy the nearest T interstitial sites along the hexagonal direction (T1T2), the second-neighbor T interstitial sites along the hexagonal direction (T1T3), the nearest T sites in the two nearest hexagonal directions (T1T4), the T and its nearest O sites (T1O), and the third-neighbor T sites along the hexagonal direction (T2T3). The corresponding formation energies of the different configurations are shown in Table 4. It is clear that the present potential reproduces the hydrogen ordering in Er, and demonstrates that the T1T3 and T1T2 configurations have the lowest and highest formation energies, respectively. The results reveal that an H pair prefers to form a configuration where the two H atoms are located at the second-neighbor T interstitial sites along the hexagonal direction, which is in agreement with the present PW91/DFT results. The properties of an ErH dimer were not considered in the fitting procedure. However, it is noteworthy that the final values of D0 = 2.85 eV and r0 = 0.18 nm are comparable to the PW91/ DFT values of 3.13 eV and 0.19 nm. Finally, the migration properties of a single H in hcp-Er were determined using the present potential in a simulation box consisting of 1440 atoms. Here, we calculate the migration energy barrier of a H atom using a drag method at a fixed volume and the atomic positions are relaxed in a hyperplane perpendicular to the vector from the initial to final positions.52 The migration energies for the diffusion paths of TT jumps perpendicular to the basal plane and TOT jumps parallel the basal plane are 0.09 and 0.52 eV, respectively, which are consistent with the PW91/DFT calculations of 0.10 and 0.53 eV for the same migration paths. B.2. ErH2. Interest in erbium hydrides was stimulated for many years because of their use as reactor moderators and also as target materials in neutron generators.2 The ErH interaction fitted in the present potential should reproduce not only the properties of H in hcp-Er, but also the bulk properties of ErH2. The results obtained with the potentials and the experimental data of ErH2 are presented in Table 5, along with the PW91/DFT calculations for comparison. As can be seen from Table 5, the lattice constant
*a
5.126
*Eb
3.679
5.129
5.122 3.68
*C11
143.4
140 ( 4
138
*C12 *C44
64.6 77.1
63 ( 2 78 ( 1
65 74
*B
90.9
Y G *EfVtet *EfFrenke Em Hoct Em Vtet
97 ( 4c
89.1 ( 1
145.3
148 ( 20c
146
139.8
58.9
60 ( 10c
58.8
55.7
89.33
1.37 0. 93
1.43 0.6b
1.21
1.20
0.70
0.62
0.56
1.39
0.98
0.92
a
Values listed are the binding energy Eb (eV/atom), the lattice parameter a (Å), the elastic constants Cij, the bulk modulus B, Young’s modulus Y and shear modulus G (GPa), the relaxed H vacancy formation energy Efvac, the formation energy EfFrenke of a Frenkel pair, and b c the migration energy Em vac of an H vacancy. Reference 2. Reference 45.
of ErH2 is consistent with the experimental data. The binding energy of ErH2 was referred to the PW91/DFT calculations, because of the lack of available experimental data. The binding energy of ErH2 is defined as Eb ¼
1 ðEN ET Þ N
ð13Þ
where N is the total number of atoms and EN and ET are the total energy of N free atoms and the total energy of a computational cell consisting of N atoms, respectively. As can be seen from Table 5, the binding energy of ErH2 calculated using the present potential is 3.68 eV/atom, which is in excellent agreement with the PW91/ DFT data (3.679 eV/atom). Also, it is clear that the present potential reproduces well the mechanical properties of bulk ErH2, such as the elastic constants, the bulk modulus, Young’s modulus and shear modulus. All H atoms are located at the tetrahedral sites in a perfect fcc ErH2. In the present work, we also consider the properties of H defects, including the formation energy of an H vacancy and a defect H pair. The H vacancy means that one H atom is taken away from the simulation box and forms one vacant tetrahedral site (Vtet). The configuration of a defect H pair consists of an H atom located at an octahedral site, with a vacancy at its originally tetrahedral site. The pair must be separated by at least one nearest neighbor distance, since an H atom is not locally stable at the octahedral site if there is a neighboring tetrahedral vacancy.9 The formation energy of an H vacancy, EfVtet, is 1.37 eV with the present potential, which agrees with the PW91/DFT value of 1.43 eV. It should be noted that the formation energy of a defect H pair (EfFrenke) predicted by the PW91/DFT method is about 1.21 eV,9 which is significantly larger than the experimental value of 0.6 eV.3 In the fitting procedure, the formation energy of a defect H pair is set to 0.90 eV, which is the mean value of the experimental result and the PW91/DFT data. The present potential gives the formation energy of a defect H pair to be 25102
dx.doi.org/10.1021/jp2090523 |J. Phys. Chem. C 2011, 115, 25097–25104
The Journal of Physical Chemistry C
ARTICLE
systems, the melting point of Er, and the diffusion barriers of H interstitials in hcp-Er. It is expected that the present potentials can be employed to investigate the dynamic behaviors in hcp-Er containing different H concentrations.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected];
[email protected].
Figure 3. Thermal expansion of the lattice parameter of a fcc ErH2 crystal, where the squares and circles represent the present results and experimental data,3 respectively.
0.86 eV, slightly smaller than the objective value of 0.90 eV. It is of interest to note that different hydrogen diffusion mechanisms appear in different H/Er stoichiometries. For H/Er > 2.0, the H atoms prefer the octahedral sites in a fcc structure, where significant hydrogen self-diffusion occurs,9 whereas, the Vtet migration is dominant for He/Er < 2.0. The migrations of octahedral interstitial H (Hoct) and Vtet may affect the stability of ErH2 during the heating process. Therefore, it is important to determine their migration energies and mechanisms. The migration barrier of Hoct (Em Hoct) is calculated to be 0.70 eV with the potential, which is close to the PW91/DFT value of 0.62.9 The migration barrier of Vtet(Em Vtet) is determined to be 1.39 eV, which is slightly larger than the PW91/DFT value of 0.98 eV.9 We found that the lower Em Vtet would lead to the collapse of the ErH2 structure. Finally the stability of ErH2 has been tested with a simulation box of 2592 atoms at zero pressure for at least 0.2 ns. Figure 3 shows the variation of the lattice parameter of ErH2 crystal as a function of temperature, where blue squares and black circles represent the values determined using the present potential and experimental results, respectively. From Figure 3, it can be seen that the change in the lattice parameter follows the behavior similar to experimental variation, i.e. the lattice parameter increases slowly at low temperatures, but shows a fast increase at higher temperatures. The present study shows that after heating the ErH2 to 700 C for 0.1 ns, the lattice change is about 0.0098 nm, which is consistent with the experimental results that suggest a very small effect on the lattice spacing after heating to 700 C for 30 min.3
IV. SUMMARY Empirical many-body interatomic potentials for ErEr and ErH interactions have been developed based on an analytical bond-order formalism, which can be used to study some important properties of an ErH system. The potentials are fitted to point defect properties, the formation energies of H interstitials and structural properties of Er and ErH systems, as obtained by experiments and our PW91/DFT calculations, using a lattice relaxation fitting approach. The present potentials can be used to describe the different phases of an ErH system with different H concentrations. In combination with Brenner’s HH potential, we demonstrate that the fitted potentials are able to reproduce the energetic and structural properties of perfect Er and ErH2 crystals, the defect properties in both Er and ErH
’ ACKNOWLEDGMENT S.M.P. and X.G.L. are grateful for the Science and Technology Foundation of China Academy of Engineering Physics (Grant No. 2009A0301015). L.Y. and X.T.Z. are grateful for the support by National Natural Science Foundation of ChinaNSAF (Grant No. 10976007) and the Fundamental Research Funds for the Central Universities (Grant No. ZYGX2009J040). F.G. is grateful for the support by the U.S. Department of Fusion Energy Science under Contract DE-AC06-76RLO 1830. ’ REFERENCES (1) Palasyuk, T.; Tkacz, M. Solid State Commun. 2005, 133, 481. (2) Grimshaw, J. A.; Spooner, F. J.; Wilson, C. G.; Mcquillan, A. D. J. Mater. Sci. 1981, 16, 2855. (3) Bonnet, J. E.; Daou, J. N. J. Appl. Phys. 1977, 48, 964. (4) Wang, Y.; Chou, M. Y. Phys. Rev. B 1994, 49, 10731. (5) Parish, C. M.; Snow, C. S.; Brewer, L. N. J. Mater. Res. 2009, 24, 1868. (6) Palasyuk, T.; Tkacz, M.; Vajda, P. Soli. State Commu. 2005, 135, 226. (7) Muller, W. M.; Blackledge, J.; Libowitz, G. G. Metal Hydrides; Academic Press: New York, 1968. (8) Palasyuk, T.; Tkacz, M. Solid State Commun. 2004, 130, 219. (9) Wixom, R. R.; Browning, J. F.; Snow, C. S.; Schultz, P. A.; Jennison, D. R. J. Appl. Phys. 2008, 103, 123708. (10) Yang, L.; Peng, S. M.; Long, X. G.; Gao, F.; Heinisch, H. L.; Kurtz, R .J.; Zu, X. T. J. Appl. Phys. 2010, 107, 054903. (11) Baskes, M. I.; Johnson, R. A. Modell. Simul. Mater. Sci. Eng. 1994, 2, 147. (12) Hu, W. Y.; Xu, H. D.; Shu, X. L.; Yuan, X. J.; Gao, B. X.; Zhang, B. W. J. Phys.: Condens. Matter 2000, 33, 711. (13) Juslin, N.; Erhart, P.; Tr€askelin, P.; Nord, J.; Henriksson, K. O. E.; Nordlund, K.; Salonen, E.; Albe, K. J. Appl. Phys. 2005, 98, 123520. (14) Tersoff, J. Phys. Rev. Lett. 1986, 56, 632. Phys. Rev. B 1988, 37, 6991. Phys. Rev. Lett. 1988, 61, 2879. Phys. Rev. B 1989, 39, 5566. (15) Brenner, D. W. Phys. Rev. B 1990, 42, 9458. (16) Pettifor, D. G.; Oleinik, I. I. Phys. Rev. B 1999, 59, 8487. (17) Wang, Z. G.; Zu, X. T.; Gao, F.; Weber, W. J. Phys. Rev. B 2008, 77, 224113. (18) Li, X. C.; Shu, X. L.; Liu, Y. N.; Gao, F.; Lu, G. H. J. Nucl. Mater. 2011, 408, 12. € € (19) BjOrkas, C.; Juslin, N.; Timko, H.; VOrtler, K.; Nordlund, K.; Henriksson, K.; Erhart, P. J. Phys.: Condens. Matter 2009, 21, 445002. (20) Albe, K.; Nordlund, K.; Averback, R. S. Phys. Rev. B 2002, 65, 195124. (21) Erhart, P.; Juslin, N.; Goy, O.; Nordlund, K.; Muller, R.; Albe, K. J. Phys.: Condens. Matter 2006, 18, 6585. (22) Mrovec, M.; Nguyen-Manh, D.; Pettifor, D. G.; Vitek, V. Phys. Rev. B 2004, 69, 094115. (23) M€uller, M.; Erhart, P.; Albe, K. J. Phys.: Condens. Matter 2007, 19, 326220. (24) Ahlgren, T.; Heinola, K.; Juslin, N; Kuronen, A. J. Appl. Phys. 2010, 107, 033516. (25) Kresse, G.; Furthmuller, J. Phys. Rev. B 1996, 54, 11169. (26) Wang, Y.; Perdew, J. P. Phys. Rev. B 1991, 44, 13298. 25103
dx.doi.org/10.1021/jp2090523 |J. Phys. Chem. C 2011, 115, 25097–25104
The Journal of Physical Chemistry C
ARTICLE
(27) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59, 1758. (28) Bl€ochl, P. E. Phys. Rev. B 1994, 50, 17953. (29) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671. (30) Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13, 5188–5192. (31) Wu, X. Z.; Wang, R.; Wang, S. F. Appl. Surf. Sci. 2010, 256, 3409. (32) Jonsson, H.; Mills, G.; Jacobsen, K. W. Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Ciccotti, G., Coker, D. F., Eds.; World Scientific: Singapore, 1998. (33) LAMMPS Molecular Dynamics Simulator, http://lammps. sandia.gov/. (34) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (35) Kittel, C. Introduction to Solid State Physics, 6th ed.; Wiley: New York, 1976; p 74. (36) Barrett, C. S.; Massalski, T. B. Structure of Metals, 3rd ed.; Pergamon: New York, 1980; p 629. (37) Brandes, E. A. Smithells Metals Reference Book, 6th ed.; Butterworths: London, 1983; p 515. (38) De Boer, F. R.; Boom, R.; Mattens, W. C. M.; Miedema, A. R.; Niessen, A. K. Cohesion in Metals: Transition Metal Alloys; North-Holland: Amsterdam, 1988. (39) http://www.webelements.com. (40) Morris, J. R.; Wang, C. Z.; Ho, K. M; Chan, C. T. Phys. Rev. B 1994, 49, 3109. (41) Balasubramanian, K.; Ma, Z. J. Phys. Chem. 1991, 95, 9794. (42) Bacon, D. J. J. Nucl. Mater. 1988, 159, 176. (43) Hu, W. Y.; Zhang, B. W.; Huang, B. Y.; Gao, F.; Bacon, D. J. J. Phys.: Condens. Matter 2001, 13, 1193. (44) Blaschko, O.; Pleschiutschning, J.; Vajda, P.; Burge, J. P.; Daou, J. N. Phys. Rev. B 1989, 40, 5344. (45) Snow, C. S.; Schultz, P.; Mattsson, T.; Bond, G.; Browning, J.; Zavadil, K.; Brumbach, M. A Hearty Hodge Podge of ErT2 Research submitted http://neutrons.ornl.gov/conf/HHIM2010/. (46) Ganchenkova, M. G.; Borodin, V. A.; Nieminen, R. M. Phys. Rev. B 2009, 79, 134101. (47) Yang, L.; Peng, S. M.; Long, X. G.; Gao, F.; Heinisch, H. L.; Kurtz, R. J.; Zu, X. T. J Phys.: Conden. Matter 2011, 23, 035701. (48) Seleskaia, T.; Osetsky, Yu. N.; Stoller, R. E.; Stocks, G. M. J. Nucl. Mater. 2007, 367370, 355. (49) Palasyuk, T.; Tkacz, M. Solid State Commun. 2004, 130, 219. (50) Blaschko, O.; Pleschiutschnig, J.; Vajda, P.; Burger, J. P.; Daou, J. N. Phys. Rev. B 1989, 40, 5344. (51) Blaschko, O.; Krexner, G.; Daou, J. N.; Vajda, P. Phys. Rev. Lett. 1985, 55, 2876. (52) Fu, C. C.; Willaime, F.; Ordejon, P. Phys. Rev. Lett. 2004, 92, 175503.
25104
dx.doi.org/10.1021/jp2090523 |J. Phys. Chem. C 2011, 115, 25097–25104
