A First-Principles Study of Hydrogen Diffusivity and Dissociation on δ

Jul 25, 2017 - The diffusivity and chemical reactivity of hydrogen on the δ-Pu (100) and (111) surfaces has been studied with density functional theo...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCC

A First-Principles Study of Hydrogen Diffusivity and Dissociation on δ‑Pu (100) and (111) Surfaces Nir Goldman* and Miguel A. Morales Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, California 94550, United States ABSTRACT: The diffusivity and chemical reactivity of hydrogen on the δ-Pu (100) and (111) surfaces has been studied with density functional theory using both spinpolarization and spin−orbit coupling calculations. Comparison of the total electronic density of states and atomic hydrogen diffusion energies indicates that spin-polarization yields accurate results, with spin−orbit coupling yielding slightly smaller barriers to hopping between adsorption sites. On the (100) surface, both sets of calculations indicate that the dissociation reaction for molecular hydrogen is highly active at ambient conditions and results in the hydrogen ions bonded within the Pu surface, similar to hydride formation. In contrast, calculations on the (111) surface indicate a lower barrier for dissociation without the formation of the hydride-like end product. Our results help quantify the effect of spin−orbit coupling on hydrogen surface chemistry on δ-Pu. In addition, we observe that the degree of resulting corrosion could depend on the specific geometry of the plutonium surface, which could have ramifications for future engineering applications.

1. INTRODUCTION Plutonium metal is highly complex from a physical and chemical perspective, with six allotropes between ambient and melting conditions.1 The metastable face-centered-cubic (fcc) δ-phase is of particular interest for engineering applications in part due to its high ductility, in contrast to the unalloyed αphase ground state, which is monoclinic and brittle.2 The δphase is stabilized by the addition of a few atomic percent of Ga3 and exhibits a number of exotic physical properties, including a negative coefficient of thermal expansion and the ability to become superconducting when alloyed.4 However, corrosion from gases in the ambient atmosphere of the application environment presents significant challenges due to the high degree of reactivity of the material and the ubiquity of these gases in experimental setups. Surface hydrogenation in particular is problematic, where the Pu-hydride product can catalyze violently exothermic oxidation reactions and even induce pyrophoricity.5 Despite these pressing issues, relatively little is known about the chemical mechanism for corrosion, and experimental studies would greatly benefit from detailed atomistic knowledge of the chemical reactivity of hydrogen on δ-Pu surfaces. Accurate modeling of the breaking and forming of bonds in condensed phases frequently requires the use of quantum simulation approaches, such as Kohn−Sham Density Functional Theory (DFT). DFT remains one of the most popular modeling methods in condensed matter physics, computational chemistry, and materials science, and has been widely used in studies of the physical and structural properties of δ-Pu.6−12 Calculations with standard DFT approaches can be inaccurate © 2017 American Chemical Society

for bulk properties such as the equilibrium volume, bulk modulus, and coefficient of thermal expansion,13 in part due to an underestimation of the correlations between the 5f electrons. However, such issues can be overcome with the inclusion of onsite Coulomb repulsion (e.g., LDA/GGA+U)9 or orbital polarization approaches,13 which can greatly improve predictions of material properties. In contrast, comparatively little effort has been made to determine surface chemistry on δ-Pu, in part due to the high computational cost from the larger system sizes needed to model the vacuum/metal surface interface. Previous studies on hydrogen adsorption and reactivity have tended to focus on atomic clusters14 or small system sizes,4,15−17 and have used the more computationally efficient spin-polarization treatment of electron spin, which explicitly neglects coupling of the electronic spin to its motion (spin−orbit coupling effects). A number of studies have additionally focused on the dissociation reaction of O2 on δ-Pu surfaces using collinear spin-polarization.18,19 Accurate treatment of electronic spin, though, is central to corrosion studies, because spin−orbit coupling effects are thought to be quite significant for plutonium, and could lower molecular dissociation barriers to exceedingly small values.14 Spin−orbit coupling calculations of atomic adsorption of carbon, nitrogen, and oxygen on the δ-Pu (111) surface showed changes of up to 0.27 eV in adsorption energies relative to spin-polarization results.20 This energy difference could Received: May 23, 2017 Revised: July 19, 2017 Published: July 25, 2017 17950

DOI: 10.1021/acs.jpcc.7b04992 J. Phys. Chem. C 2017, 121, 17950−17957

Article

The Journal of Physical Chemistry C

values have been tuned to yield accurate bulk moduli, albeit with relatively large energy convergence criteria (10−3 eV).32 However, our goal is to create an accurate minimum energy geometry for different δ-Pu slabs for the study of surface chemistry. In this respect, hydrogen dissociation chemistry did not show any measurable dependence on the exact U − J values used (discussed below). Care was taken to minimize each geometry to as close to an antiferromagnetic surface (AFM) as possible.7,8,33 The magnetism of δ-Pu is complicated,34 and has been the subject of a number of recent experimental and theoretical studies.10,11,35 Regardless, an AFM state yields structural properties comparable to experiments within the level of theory employed in this work (e.g., GGA+U), and is lower in energy than nonAFM states.3 We note that there are any number of spin arrangements that will yield the AFM material, although these structures tend to be nearly iso-energetic.6,13 We have confirmed this by computing the AFM total energy for an optimized 32 atom δ-Pu lattice for one geometry where the magnetic spins have layers in the ⟨100⟩ direction of alternating positive and negative values, and a second where the Pu atoms within a given layer have alternating sign. In this case, we observe an energy difference of 9.42 meV/atom (0.217 kcal/ mol) between each spin configuration for the SP result, with the geometry with antiparallel spins within the Pu layer having the slightly lower energy. Similar results were found for the SP computed H2 dissociation energy map on the (100) surface, where we observed a mean absolute deviation of 0.05 eV between energies for specific adsorption configurations, with deviations of only 0.01 eV for specific reaction barrier heights. In addition, the spin geometry with alternating planes of magnetic moments tended to yield a nonzero net magnetic moment for our H2 dissociation calculations on the (100) and (111) surfaces. As a result, in this work we report surface chemistry results for the spin geometry with alternating signs within each layer only.

prove significant for atomic diffusion rates at room temperature (298 K ≈ 0.026 eV). Inclusion of noncollinear spin−orbit coupling has also been found to be significant for bulk δplutonium hydriding and dehydriding energies.21 These findings are in contrast to studies on the α-Pu phase, where spin−orbit coupling had a negligible effect on hydrogen atomic adsorption.22 Nonetheless, the importance of spin−orbit coupling in determination of the chemical reactivity of atmospheric gases on plutonium surfaces has not been fully assessed, and its effects relative to spin-polarization calculations remain largely unquantified. To this end, we have computed the surface adsorption of atomic and molecular hydrogen and the energetics of the dissociation reaction using both collinear spin-polarization and noncollinear spin−orbit coupling. In addition to its significance regarding plutonium corrosion, hydrogen adsorption and dissociation are the simplest prototypical chemical reactions that can yield fundamental insight into the complexity of any plutonium surface chemical reactivity. We have computed the electronic density of states for bulk δ-Pu as well as the potential diffusivity of atomic hydrogen on the (100) and (111) surfaces as a way of directly comparing results from both methods. We have then computed the manifold of H2 dissociation pathways on both surfaces and have determined both reaction mechanisms as well as the resulting energetic barriers. Our results quantify the differences between spin-polarization and spin−orbit coupling for these systems, and provide a more detailed prediction for hydrogen diffusion and dissociation reaction paths than previously available. These findings could have significant ramifications for future experiments and calculations on corrosive reactions on plutonium-based materials.

2. COMPUTATIONAL DETAILS All DFT calculations were performed with the VASP code.23−25 We label collinear spin-polarization calculations as SP and noncollinear spin−orbit coupling calculations as SOC for the sake of brevity. For both SP and SOC approaches, calculations were performed with projector-augmented wave function (PAW) pseudopotentials26,27 and the Perdew, Burke, and Ernzerhof exchange correlation functional (PBE).28 Partial occupancies of the electronic states were set with a fourth-order Methfessel−Paxton smearing29 with a width of 0.3 eV. All energies reported here are without the electronic entropic effects. Wave function convergence criteria of 10−4 eV were found to be sufficient for the SP calculations, although SOC calculations required a more stringent criteria of 10−5 eV due to the presence of numerable local minima. Calculations on an optimized δ-Pu four atom unit cell with a range of cutoffs up to 600 eV indicated that our results are sufficiently converged at 318 eV. For this work, we have used a cutoff of 400 eV for all calculations. Additional calculations with a 600 eV cutoff for H2 adsorption on an optimized δ-Pu (100) surface agree within 0.1 eV. To accurately account for the effective correlation between forbitals, we have used the GGA+U approach of Dudarev et al.,30 as implemented in VASP. The value of the effective onsite Hubbard correction U − J was varied to reproduce a cell volume as close as possible to the experimental result.31 For the SP calculations, this resulted in a U − J value of 1.30 eV (U = 2.05 eV and J = 0.75 eV). For the SOC calculations, we used a U − J value of 1.00 eV (U = 1.75 eV and J = 0.75 eV). These value are somewhat smaller than previous results, where U − J

3. RESULTS AND DISCUSSION To compare bulk properties from SP and SOC methods, we first performed a lattice optimization on a unit cell of four atoms with a Monkhorst−Pack mesh36 of 10 × 10 × 10, although the lattice vectors were already converged within 0.1 Å with a mesh as small as 6 × 6 × 6 for both SP and SOC calculations. Forces were converged to 0.01−0.05 eV/Å. The resulting lattice parameters (Table 1) indicate a slight distortion Table 1. Lattice Parameters for δ-Pu from SP and SOC Optimizations method 31

expt. SP SOC

a (Å)

b (Å)

c (Å)

volume (Å3/atom)

4.637 4.581 4.619

4.808 4.677

4.581 4.619

24.93 25.23 24.95

of the initial fcc lattice into a tetragonal lattice, similar to previous results,6,12 with the SOC result slightly closer to experiment.31 A number of fcc actinide compounds exhibit tetragonal lattice distortion induced by antiferromagnetism.37 A cubic phase for δ-Pu has been observed in low temperature Xray diffraction experiments,38 although the tetrahedral phase is known to be stabilized at 100 K.7 The optimized δ-Pu lattice was then replicated in each direction to yield a simulation cell of 32 Pu atoms, or four 17951

DOI: 10.1021/acs.jpcc.7b04992 J. Phys. Chem. C 2017, 121, 17950−17957

Article

The Journal of Physical Chemistry C layers of eight atoms each, and was reoptimized using a k-point mesh of 5 × 5 × 5 for both SP and SOC calculations. We then computed the total electronic density of states (DOS) from both methods for the larger optimized simulation cell. Here, we have normalized the DOS to a value of one to facilitate comparison (Figure 1). We observe that the SOC and SP

states (PDOS) for both SP and SOC methods for each surface indicates similar features (Figure 2), where the vast majority of the energy states around the Fermi energy are due to f-orbital interactions. We have computed the average f-electron band energy for the (100) surface to be −0.52 eV for the SP calculation and −0.61 eV for SOC, and values for the (111) surface of −0.26 for the SP calculation and −0.57 for SOC. These results indicate the likelihood of similar adsorption energies and surface chemistry between the two methods, given their relatively similar values as well as the fact that all are only slightly below the Fermi energy.39 Inspection of the (100) and (111) surfaces (Figure 3) shows distinct differences in surface geometries that ultimately lead to significant differences in the resulting chemical properties. The (100) surface contains relatively small interstitial sites between surface Pu atoms, and exhibits a clear fcc-like ABAB stacking, where atoms residing in the first subsurface layer are centered essentially directly below the center of the interstitial opening. In contrast, the (111) surface contains much larger interstitial openings, with the layers now adopting a more hexagonal ABCABC stacking, and the first subsurface layer is aligned far from the center. We have computed energy maps for possible hydrogen atomic surface diffusion by performing energy grid calculations, where a single hydrogen atom was placed 2.0 Å above the surface over a grid of points, and then had its height optimized relative to the surface while its two transverse coordinates were held fixed. Because of the large mass discrepancy with hydrogen, the plutonium atoms were frozen in their optimized slab positions. In doing so, we have created a grid of all possible adsorption energies and subsequent diffusion pathways on each crystalline surface. Here, grid points were placed in increments of approximately 0.5 Å along each transverse axis, spanning one lattice displacement (e.g., ∼4.5 Å), for a total of 16 calculations. On each surface and for each calculation method, optimal configurations were confirmed by computing single point energies of H atoms displaced by ±0.1−0.2 Å along the surface normal. For the (100) surface, we observe very similar results between SP and SOC calculations (Figure 4). In this work, we have defined the adsorption energy as Eads = Eslab+adsorbate − (Eslab + Eadsorbate). Both SP and SOC results indicate qualitatively the same surface mobility pathways, with the hydrogens exhibiting preferential adsorption above the interstitial opening (center sites), and higher energies for adsorption directly above surface Pu atoms. Specific adsorption energies and approximate distances from the surface for specific positions are shown in Table 3. For the center sites, the SP calculations indicate an adsorption energy of −3.00 eV, and SOC energies of −2.81 to −2.91 eV. In this case, the hydrogen atom adsorbed directly into the surface layer of the Pu lattice, corresponding to the first step in forming a plutonium hydride, where hydrogen embeds itself directly into the plutonium surface layer. The hydrogen experiences slightly less favorable adsorption directly between two neighboring surface Pu atoms (bridging sites), with the SP calculations yielding values of −2.59 to −2.65 eV, and SOC with values from −2.48 to −2.68 eV. Here, the hydrogen atom remained ∼1.5 Å above the surface for the SP calculations and ∼1.6 Å above for the SOC calculations. For adsorption directly above the Pu atoms (top sites), the hydrogen atom exhibited energies of −1.87 to −1.89 eV for the SP results, and energies of −1.76 to 1.91 eV for the SOC calculations. For comparison, previous results for the

Figure 1. Total electronic density of states computed for δ-Pu using spin−orbit coupling and spin-polarization (spin-up contributions only). In each case, the density of states has been recentered so the Fermi energy equals zero.

results are generally fairly similar. The SOC result exhibits a higher density of states at about −2 eV, which is compensated for by a slightly lower density of states at approximately −0.5 eV. However, both methods yield similar densities of state around the Fermi energy (centered at 0 eV). The (100) surface slab was then created by placing a 30 Å vacuum above the system in the ⟨100⟩ direction, and optimizing the atomic positions of the top two layers while freezing the bottom two. A k-point mesh of 1 × 5 × 5 was found to be sufficient for the SP calculations, with the smaller mesh value lying in the direction of the surface. For the sake of computational efficiency, we reduced the k-point mesh for SOC calculations to 1 × 4 × 4. To create an orthorhombic supercell, an orthorhombic unit cell of six atoms with the (111) face aligned normal to one of the principle Cartesian axes was first optimized using the same (U − J) values as for the (100) surface and a k-point mesh of 10 × 10 × 10. The cell was then replicated in the orthonormal directions to yield a simulation cell of 24 Pu atoms (6 atomic layers), with the vacuum layer placed above the (111) surface. Again, 1 × 5 × 5 and 1 × 4 × 4 k-point meshes were found to be sufficient for the SP and SOC calculations, respectively. Supercell dimensions for all surfaces and methods are listed in Table 2. The projected density of Table 2. Supercell Dimensions for SP and SOC Optimizations for Both (100) and (111) Surfacesa (100) surface (111) surface

method

lx (Å)

ly (Å)

lz (Å)

SP SOC SP SOC

39.1610 39.1570 5.8947 5.7593

9.6166 9.7226 46.2318 46.0758

9.1610 9.1570 6.3789 6.4405

a

For the (100) system, the surface was aligned normal to the lx supercell vector, whereas for the (111) system, the surface was aligned normal to ly. 17952

DOI: 10.1021/acs.jpcc.7b04992 J. Phys. Chem. C 2017, 121, 17950−17957

Article

The Journal of Physical Chemistry C

Figure 2. Projected density of states for SP and SOC calculations. In each case, f-orbital contributions are colored red, d-orbital contributions are orange, p-orbital contributions are blue, and s-orbital contributions are black, with the total DOS shown as a thin black line for reference. For the SP case, the spin-up contribution is shown as positive and the spin-down as negative.

corresponding energy from the SP result. Regardless, the range of both sets of results covers about 0.9 eV, with the SP calculations showing adsorption energies between ∼ −1.7 and −2.6 eV, and the SOC result between ∼ −1.9 and −2.8 eV. In this case, both SP and SOC calculations show the highest hydrogen atom adsorption energies directly above the surface Pu atoms, similar to the (100) surface. However, the hydrogen atoms appear to be relatively mobile on this surface, with adsorption energies that vary by only approximately 0.3 eV, excluding the surface Pu atoms, where adsorption is least favorable. There is some variation between the SP and SOC results for adsorption in the center sites on the surface. However, both sets of calculations show preferential adsorption for configurations where the H atom is located equidistant from three Pu surface atoms, about 1.5 Å above the surface. We compute nearly identical adsorption energies and surface distances in hcp-type (located directly above the first subsurface Pu atom) versus fcc-type (located above the second subsurface atom) center sites for both methods. Previous calculations indicate similar relative energies between adsorption at different surface sites with similar H atomic surface distances.15 However, the adsorption energies computed there are up to 1 eV lower in energy, possibly due to the previously mentioned small system sizes and unoptimized surface geometry. Given the relative similarities between the SP and SOC results, we have proceeded to compute H2 dissociation on both crystalline surfaces with the SP methodology, primarily. A full relaxation (surface normal and transverse directions) of the H2 molecule and top two layers of the Pu surface yielded

Figure 3. Renderings of the (100) and (111) surfaces of δ-Pu surfaces computed from noncollinear calculations using spin−orbit coupling. Surface Pu atoms are completely opaque, the first subsurface layer is colored red, and subsequent subsurface atoms are semitransparent. The (100) surface has significantly smaller interstitial openings that contribute to the H2 dissociation chemistry.

(100) surface15 exhibited similar distances for adsorption in the bridging sites and top sites above Pu atoms. However, adsorption in the center site was shown to take place 1.13 Å above the surface, in contrast to the surface embedding of hydrogen seen in this work. In addition, all of the adsorption energies differed from the results shown here by up to 0.5 eV. These discrepancies are likely due in part to the significantly smaller simulation cell of previous work (six atoms only), which also did not use a slab-optimized surface geometry. We observe somewhat different surface mobility properties on the (111) surface for both sets of calculations (Figure 5). We note that the values of the adsorption energies for the SP and SOC are shifted slightly with respect to each other with the SOC results yielding values that are ∼0.2 eV lower than the 17953

DOI: 10.1021/acs.jpcc.7b04992 J. Phys. Chem. C 2017, 121, 17950−17957

Article

The Journal of Physical Chemistry C

Table 3. Selected Hydrogen Atomic Adsorption Energies and Approximate Distance from the Surface, As Compared to Spin Polarized Results from Reference 15a surface (100)

site top

bridge

center

(111)

top

bridge

hcp-type

fcc-type

method

energy (eV)

distance (Å)

SP SOC ref 15 SP SOC ref 15 SP SOC ref 15 SP SOC ref 15 SP SOC ref 15 SP SOC ref 15 SP SOC ref 15

−1.873 −1.913 −2.122 −2.656 −2.684 −3.141 −3.016 −2.912 −3.467 −1.745 −2.022 −2.517 −2.457 −2.547 −3.126 −2.509 −2.606 −3.450 −2.487 −2.596

2.14 2.04 2.03 1.49 1.36 1.56 0.00 0.00 1.13 2.03 2.09 2.01 1.57 1.52 1.57 1.53 1.52 1.42 1.49 1.48

“Top” corresponds to adsorption directly above a surface Pu atom, “bridge” to directly between two surface Pu atoms, and “center” to positions centered between four surface Pu atoms on the (100) surface. Results from ref 15 do not include adsorption in the fcc-type site for the (111) surface.

a

Figure 4. Maps of hydrogen atomic diffusion energies on the δ-Pu (100) surface. Plutonium surface atoms (labeled) are located in areas of less favorable adsorption. For each set of calculations, the H atom experiences preferential adsorption (dark blue regions) in the interstitial openings. The apparent slight breaking of symmetry in the energy surfaces is likely due to the presence of nearly iso-energetic minima on the AFM surface, where each geometry optimization calculation has resulted in slightly different spin geometries, which can yield slight differences in adsorption energies.

bottom H atom was kept positioned within the (100) center site, and the top H atom was translated along a grid of points spaced by 0.3−1.0 Å in each transverse direction. This resulted in over 50 optimization calculations, which could be performed efficiently in parallel. We have tested the validity of our constrained dynamics pathway by performing calculations where both hydrogens were placed 2.0 Å above the surface with a dissociative bond length (e.g., 1.5 Å), with all of their cartesian coordinates optimized. In this case, we observe an identical mechanism where the hydrogen above the center site embeds itself within the Pu surface without any intermediate energy minimum configurations. This is consistent with screening effects from the Pu surface, which results in relatively short-ranged interactions between hydrogen atoms. Nonetheless, our approximate reaction surface does not take into account more complex molecular displacements, and hence represents a possible upper limit to the actual reaction barrier. A more accurate determination of minimum reaction barriers would involve full inspection of the six-dimensional potential energy surface, which is beyond the scope of the current work. The dissociation minimum energy pathway could have been computed more precisely through use of the nudged elastic band method (NEB),40,41 which is the possible subject of future work. Similarly, saddle-point search approaches such as the DIMER method42 and a vibrational frequency analysis can be used to precisely identify transition states. However, these methods can require extremely high convergence of the system energies to yield accurate atomic forces. This problem is compounded by the fact that the H2 reaction barrier is relatively low and of magnitude similar to that of other accessible reaction modes like molecular diffusion. In contrast, the energy

preferential adsorption in the center site with an Eads value of −0.04 eV. We then determined the H2 diffusion energy maps on both the (100) and the (111) surfaces in a similar fashion where the H2 geometry was optimized relative to the surface normal over a grid of points in the transverse directions, spaced by ∼1.0 Å over a range of 4.5 Å (25 points, total). For the (100) surface, we found that the unreacted H2 molecule generally optimizes to a vertical position, where the center of mass is approximately 3.6 Å above the surface. An optimization calculation with the hydrogen molecule initially in a horizontal orientation resulted in the H2 molecule oriented 19.8° from the vertical axis, which was only 0.011 eV higher in energy than the completely vertical configuration. Because of the variability resulting from different final optimized AFM spin geometries, we find it instructive to compare relative energies for adsorption for a single Pu spin arrangement rather than absolute values of adsorption energy. We observe that the H2 molecule has potentially high surface mobility, with molecular adsorption energies that vary over a range of 0.03 eV, similar to previous results.4,16 We mapped out the dissociation potential energy surface by starting with the H2 in an optimized geometry above the center site. We then performed surface normal optimizations where the positions of both H atoms were held fixed in the transverse directions, and only the position of the H atom relative to the Pu surface was allowed to vary. For these calculations, the 17954

DOI: 10.1021/acs.jpcc.7b04992 J. Phys. Chem. C 2017, 121, 17950−17957

Article

The Journal of Physical Chemistry C

Figure 6. H2 dissociation energy pathways on the (100) surface. The top panel shows the energy map for the dissociation reaction, and the bottom panel shows the height of the second H atom relative to the plutonium surface. The dashed line is used as a guide for our computed dissociation pathway. The dissociation reaction has an approximate energetic barrier of 0.22 eV, with the formation of a hydride-like end state.

Figure 5. Maps of hydrogen atomic diffusion energies on the δ-Pu (111) surface. Again, plutonium surface atoms are located in areas of less favorable adsorption. On this surface, the H atom experiences preferential adsorption (dark blue regions) when surrounded by three surface Pu atoms and large energetic barriers for hydride formation. Again, the apparent slight breaking of symmetry in the energy surfaces is likely due to the presence of nearly iso-energetic minima on the AFM surface.

bisects the ⟨010⟩ and ⟨001⟩ axes to land in the global energetic minimum configuration embedded within the Pu surface. We have validated our use of the SP approach for dissociation chemistry by computing SOC single point energies for 30 of the SP geometries determined for the (100) surface. We find that this yields total energies that have a mean absolute deviation of only 0.05 eV from the SP result, with an average surface normal force on the H atoms of less than 0.05 eV/Å, i.e., within our convergence criteria. Our SOC calculations also yield reaction energies nearly identical to those from the SP method, with a barrier value for the minimum energy path of 0.22 eV, a barrier value for dissociation along the ⟨010⟩ axis of 0.76 eV, and a total reaction energy of −1.10 eV. We have also tested the effect of the (U − J) value on the SP computed dissociation chemistry by reoptimizing the same geometries but with varying the value of (U − J) from 0 to 4 eV. We find that the results are essentially invariant to this parameter, where optimized geometries vary by only several hundredths of an angstroms and energies by several thousandths of an electronvolt (∼0.1 kcal/mol). Finally, we have computed the potential diffusion and dissociation pathways for H2 on the (111) surface (Figure 7). Geometry optimization yielded preferential adsorption in the hcp-type site with an energy of −0.16 eV. Similar to the (100) result, we find that H2 molecular adsorption is essentially isoenergetic on this surface, varying by roughly 0.04 eV over our energy grid (spanning one lattice spacing in the transverse directions, for a total of 30 grid points). This again suggests a high degree of surface mobility under ambient conditions. The center of mass distance of the molecule from the surface when adsorbed in the hcp-type center site varies from 3.2 to 3.3 Å,

grid approach used here is very general and computationally efficient, and yields reasonable upper limits for the reaction barriers and pathways. Our results indicate a complex manifold of reaction pathways (Figure 6). We observe a dissociation barrier of ∼0.22 eV for the approximate minimum energy path with a total energy of reaction of approximately −1.22 eV. These results are consistent with previous predictions of the Pu hydriding process being entirely exothermic17,43 and diffusion limited44 under ambient conditions. At the apparent reaction barrier, we observe a molecular axis that is parallel to the crystal surface. Our computed reaction barrier is significantly lower than the previously published result of 0.778 eV,16 which is likely due to the small simulation cell and possibly insufficient sampling of phase space in that study. We note that we compute a dissociation energy of approximately 0.80 eV when the reaction occurs directly along either the ⟨010⟩ or the ⟨001⟩ axes. The bottom hydrogen embeds itself directly into the Pu surface once the dissociation barrier has been crossed by the top atom, regardless of the reaction pathway. For the minimum energy path, the dissociating (originally top-most) H atom will then move along the ⟨010⟩ axis to a first minimum, approximately 1.5 Å above the Pu surface in the bridging site and approximately −0.85 eV lower in energy than the starting configuration. It will then have to overcome a dissociation barrier of roughly 0.20 eV to diffuse along a direction that 17955

DOI: 10.1021/acs.jpcc.7b04992 J. Phys. Chem. C 2017, 121, 17950−17957

Article

The Journal of Physical Chemistry C

ensure that the system remains antiferromagnetic to achieve the lowest energy configurations. However, the specific spin arrangement as well the exact values of the GGA+U correction are relatively unimportant in terms of dissociation energies. We find that both methods yield fairly similar results, with spin− orbit coupling enhancing adsorption by up to several tenths of an electronvolt, but having smaller effects on the dissociation energy barrier. Dissociation for H2 molecules adsorbed on the (100) surface is likely rapid under ambient temperature, where the final state yields hydrogen atoms that have embedded into the Pu surface. This represents the first step in Pu corrosion and the initial stage of formation of plutonium hydride. In contrast, the (111) surface exhibits slightly lower dissociation barriers but high barriers to surface embedding. Hence, this surface appears to be somewhat corrosion resistant in the presence of hydrogen. Overall, our results indicate a straightforward method for computing Pu surface chemistry, where the system can be optimized with the more computationally efficient spin-polarization method and further refined with spin−orbit coupling as needed. Our efforts can be extended to a more complete mapping of plutonium corrosion in the presence of additional atmospheric gases, which could yield advances in the way Pu-containing devices are machined for future applications.

■

Figure 7. H2 dissociation energy pathways on the (111) surface, with top and bottom panels indicating the same as Figure 6. Again, the dashed line is used as a guide for our computed dissociation pathway. The dissociation reaction has an approximate energetic barrier of 0.12 eV, with both H atoms remaining ∼1.3 Å above the surface.

AUTHOR INFORMATION

Corresponding Author

*Phone: (925) 422-3994. E-mail: [email protected]. ORCID

with the molecular axis at an angle of ∼31.1° relative to the Pu surface. We then started with the H2 molecule over the optimal hcp-type center site and computed the dissociation energy maps, as before. Here, we have reduced the grid spacing to ∼0.25 Å in the transverse directions, resulting in a total of 160 surface calculations. We compute a dissociation barrier of 0.12 eV for the minimum energy path, which is lower than the previously reported result of 0.305 eV.4 The molecular geometry at the approximate reaction maximum again yielded a hydrogen molecular axis parallel to the crystal surface. Calculation of the dissociation energy barrier with the SOC method with geometries determined from SP yielded a value of 0.09 eV. For the minimum energy path, the bottom H atom settles into a configuration in the hcp-type center site, approximately 1.3 Å above the surface once the top atom has crossed the approximate reaction barrier. The top H atom will first settle into an energetic minimum that is approximately −0.74 eV lower in energy than the system starting point at a neighboring fcc-type site. It will then cross a barrier of approximately 0.25 eV before settling into a global energetic minimum configuration at a neighboring hcp-type, which yields a total energy of reaction of approximately −0.94 eV. Although the energetic barrier for dissociation is slightly lower than our result for the (100) surface, the H atoms are still unable to penetrate the Pu surface itself. Subsequent to dissociation, each H atom experiences low barriers to diffusion, similar to the atomic energy grid result.

Nir Goldman: 0000-0003-3052-2128 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. We thank Bill Mclean, Kerri Blobaum, Pat Allen, and Albert Loui for helpful discussions.

■

REFERENCES

(1) Atta-Fynn, R.; Ray, A. K. Does hybrid density functional theory predict a non-magnetic ground state for δ-Pu. EPL 2009, 85, 27008. (2) Arsenlis, A.; Wolfer, W. G.; Schwartz, A. J. Change in flow stress and ductility of δ-phase Pu-Ga alloys due to self-irradiation damage. J. Nucl. Mater. 2005, 336, 31−39. (3) Hernandez, S. C.; Schwartz, D. S.; Taylor, C. D.; Ray, A. K. Ab initio study of gallium stabilized δ-plutonium alloys and hydrogenvacancy complexes. J. Phys.: Condens. Matter 2014, 26, 235502. (4) Huda, M. N.; Ray, A. K. Molecular hydrogen adsorption and dissociation on the plutonium (111) surface. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 085101. (5) Sun, B.; Liu, H.; Song, H.; Zhang, G.; Zheng, H.; Zhao, Z.-G.; Zhang, P. The different roles of Pu-oxide overlayers in the hydrogenation of Pu-metal: An ab initio molecular dynamics study based on van der Waals density functional (vdW-DF)+U. J. Chem. Phys. 2014, 140, 164709. (6) Sö derlind, P.; Landa, A.; Sadigh, B. Density-functional investigation of magnetism in δ-Pu. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 66, 205109. (7) Söderlind, P.; Sadigh, B. Density-Functional calculations of α, β, γ, δ, δ′, and ϵ Plutonium. Phys. Rev. Lett. 2004, 92, 185702.

4. CONCLUSIONS We have performed extensive calculations of hydrogen diffusion and dissociation on the (100) and (111) surface of δ-Pu using both collinear spin-polarization and noncollinear spin−orbit coupling methods. We find that great care must be taken to 17956

DOI: 10.1021/acs.jpcc.7b04992 J. Phys. Chem. C 2017, 121, 17950−17957

Article

The Journal of Physical Chemistry C

(32) Chen, S. P. Correlation-induced anomalies and extreme sensitivity in fcc Pu. Philos. Mag. 2009, 89, 1813−1822. (33) Wang, Y.; Sun, Y. First-principles thermodynamic calculations for δ-Pu and ϵ-Pu. J. Phys.: Condens. Matter 2000, 12, L311−316. (34) Lashley, J. C.; Lawson, A.; McQueeney, R. J.; Lander, G. H. Absence of magnetic moments in plutonium. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 054416. (35) Janoschek, M.; Das, P.; Chakrabarti, B.; Abernathy, D. L.; Lunsden, M. D.; Lawrence, J. M.; Thompson, J. D.; Lander, G. H.; Mitchell, J. N.; Richmond, S.; et al. The valence-fluctuating ground state of plutonium. Sci. Adv. 2015, 1, e1500188. (36) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188−5192. (37) Lander, G.; Mueller, M. Magnetically induced lattice distortions in actinide compounds. Phys. Rev. B 1974, 10, 1994. (38) Jacquemin, J.; Lallement, R. Plutonium and Other Actinides: Proceedings of the Fourth International Conference on Plutonium and Other Actinides; 1970; pp 616−622. (39) Hammer, B.; Norskov, J. K. Why gold is the noblest of all the metals. Nature 1995, 376, 238−240. (40) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 2000, 113, 9901. (41) Goldman, N.; Browning, N. D. Gold Cluster Diffusion Kinetics on Stoichiometric and Reduced Surfaces of Rutile TiO2(110). J. Phys. Chem. C 2011, 115, 11611−11617. (42) Henkelman, G.; Jónsson, H. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives. J. Chem. Phys. 1999, 111, 7010. (43) Taylor, C. D.; Francis, M. F.; Schwartz, D. S. Thermodynamics, structure, and charge state of hydrogen-vacancy complexes in δplutonium. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 214114. (44) Stakebake, J. L. Kinetics for the reaction of hydrogen with plutonium powder. J. Alloys Compd. 1992, 187, 271−283.

(8) Söderlind, P. Quantifying the importance of orbital over spin correlations in δ-Pu within density-functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 085101. (9) Liu, T.; Cai, T.; Gao, T.; Li, G. The electronic and structural properties of δ-Pu and PuO from the LSDA (GGA)+U method. Phys. B 2010, 405, 3717−3721. (10) Söderlind, P.; Zhou, F.; Landa, A.; Klepeis, J. E. Phonon and magnetic structure in δ-plutonium from density functional theory. Sci. Rep. 2015, 5, 15958. (11) Migliori, A.; Söderlind, P.; Landa, A.; Freibert, F. J.; Maiorov, B.; Ramshaw, B. J.; Betts, J. B. Origin of the multiple configurations that drive the response of δ-plutonium’s elastic moduli to temperature. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 11158−11161. (12) Hernandez, S.; Freibert, F.; Wills, J. Density functional theory study of defects in unalloyed δ-Pu. Scr. Mater. 2017, 134, 57−60. (13) Söderlind, P. Ambient pressure phase diagram of plutonium: A unified theory for α-Pu and δ-Pu. Europhys. Lett. 2001, 55, 525−531. (14) Balasubramanian, K.; Felter, T. E.; Anklam, T. A.; Trelenberg, T. W.; W, M., II Atomistic level relativistic quantum modelling of plutonium hydrogen reaction. J. Alloys Compd. 2007, 444−445, 447452. (15) Huda, M. N.; Ray, A. K. A density functional study of atomic hydrogen adsorption on plutonium layers. Phys. B 2004, 352, 5−17. (16) Huda, M. N.; Ray, A. K. An ab initio study of H2 interaction with the Pu (100) surface. Phys. B 2005, 366, 95−109. (17) Taylor, C. D.; Hernandez, S. C.; Francis, M. F.; Schwartz, D. S.; Ray, A. K. Hydrogen trapping in δ-Pu: insights from electronic structure calculations. J. Phys.: Condens. Matter 2013, 25, 265001. (18) Huda, M. N.; Ray, A. K. A density functional study of molecular oxygen adsorption and reaction barrier on Pu (100) surface. Eur. Phys. J. B 2005, 43, 131−141. (19) Huda, M. N.; Ray, A. K. Ab Initio Study of Molecular Oxygen Adsorption on Pu (111) Surface. Int. J. Quantum Chem. 2005, 105, 280−291. (20) Atta-Fynn, R.; Ray, A. K. Ab initio full-potential fully relativistic study of atomic carbon, nitrogen, and oxygen chemisorption on the (111) surface of δ-Pu. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 195112. (21) Yang, Y.; Zhang, P. Hydriding and dehydriding energies of PuHx from ab initio calculations. Phys. Lett. A 2015, 379, 1649−1653. (22) Islam, M. F.; Ray, A. K. An ab initio study of hydrogen adsorption on the (020) surface of α-Pu. Phys. Status Solidi B 2011, 248, 193−202. (23) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561. (24) Kresse, G.; Hafner, J. Ab initio molecular dynamics simulation of the liquid-metal-amorphous-semiconductor transition in germanium. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 14251−14271. (25) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (26) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (27) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758−1775. (28) Perdew, J. P.; Burke, K.; Enzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (29) Methfessel, M.; Paxton, A. T. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 40, 3616−3621. (30) Dudarev, S.; Botton, G.; Savrasov, S.; Humphreys, C.; Sutton, A. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 57, 1505. (31) Wallace, D. C. Electronic and phonon properties of six crystalline phases of Pu metal. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 15433−15439.

■

NOTE ADDED AFTER ASAP PUBLICATION This paper was published on August 8, 2017. Figure 2 has been corrected and the paper was re-posted on August 9, 2017.

17957

DOI: 10.1021/acs.jpcc.7b04992 J. Phys. Chem. C 2017, 121, 17950−17957