Next Article in Journal
Transfer Learning-Assisted Evolutionary Dynamic Optimisation for Dynamic Human-Robot Collaborative Disassembly Line Balancing
Next Article in Special Issue
Unified Analytic Melt-Shear Model in the Limit of Quantum Melting
Previous Article in Journal
Special Issue “Application of Non-Linear Dynamics”
Previous Article in Special Issue
Ab Initio Phase Diagram of Chromium to 2.5 TPa
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Estimates of Quantum Tunneling Effects for Hydrogen Diffusion in PuO2

1
Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
2
Department of Chemical Engineering, University of California, Davis, CA 95616, USA
3
Department of Chemical Engineering, University of Michigan, Ann Arbor, MI 48109, USA
*
Author to whom correspondence should be addressed.
Current address: Department of Chemical and Biomolecular Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA.
Appl. Sci. 2022, 12(21), 11005; https://doi.org/10.3390/app122111005
Submission received: 6 October 2022 / Revised: 25 October 2022 / Accepted: 26 October 2022 / Published: 30 October 2022
(This article belongs to the Special Issue Feature Paper Collection in Section Materials)

Abstract

:
We detail the estimation of activation energies and quantum nuclear vibrational tunneling effects for hydrogen diffusion in PuO 2 based on Density Functional Theory calculations and a quantum double well approximation. We find that results are relatively insensitive to choice of exchange correlation functional. In addition, the representation of spin in the system and use of an extended Hubbard U correction has only a small effect on hydrogen point defect formation energies when the PuO 2 lattice is held fixed at the experimental density. We then compute approximate activation energies for transitions between hydrogen interstitial sites seeded by a semi-empirical quantum model and determine the quantum tunneling enhancement relative to classical kinetic rates. Our model indicates that diffusion rates in H/PuO 2 systems could be enhanced by more than one order of magnitude at ambient conditions and that these effects persist at high temperature. The method we propose here can be used as a fast screening tool for assessing possible quantum nuclear vibrational effects in any number of condensed phase materials and surfaces, where hydrogen hopping tends to follow well defined minimum energy pathways.

1. Introduction

Plutonium oxide (PuO 2 ) is a stable ceramic that forms spontaneously from the metal in the presence of atmospheric oxygen or water [1]. PuO 2 has a number of industrial uses, including as a component in mixed oxide fuels for nuclear power reactors [2] and as a fuel source for deep-space spacecrafts [3]. Similar to many transition metals, the oxide layer on plutonium-based materials acts as a protection against corrosion due to atmospheric gases [4]. In particular, PuO 2 is thought to help prevent nucleation of hydriding due to the poor solubility of hydrogen in the lattice [5]. Ab initio molecular dynamics studies have observed that hydrogen diffusion could be slowed by the formation of O–H bonds within the bulk [6]. Despite this, plutonium metal is still highly sensitive to the presence of adsorbed hydrogen, where the Pu-hydride product can form spontaneously [7,8] and strongly exothermically [9,10,11]. Hydriding studies would benefit from atomic-scale calculations to help elucidate the initial stages of this process, including hydrogen absorption and diffusion in the bulk.
In this regard, quantum mechanical calculations using Density Functional Theory (DFT) hold promise as a method to help determine many observables for the H/PuO 2 system. DFT has been shown to yield accurate results for many bulk PuO 2 physical and mechanical properties, including the bulk modulus and band gap [12,13]. However, many challenges with performing these calculations remain. To the best of our knowledge, there are no systematic studies to determine the best choice of exchange-correlation functional for hydrogen-containing PuO 2 systems. In addition, use of an onsite Hubbard-like correction (e.g., GGA + U) [14] combined with calculation of spin effects [15] is required to yield an accurate lattice constant and band gap. Collinear spin polarization calculations tend to be more computationally efficient [16], but neglect spin–orbit coupling and relativistic effects that are known to be significant for plutonium [13]. Both GGA + U and spin–orbit coupling are generally computationally cumbersome due to their creation of numerous local minima in the electron density [17,18,19], which make it difficult to find the ground-state of the system. The magnetic ground-state of PuO 2 is itself somewhat unknown, with both experiments and DFT calculations yielding a wide range of possible spin states [12,13], and the search for the true ground-state could require placing constraints on the orbital occupations [18,19]. In contrast, the 1s orbitals from the hydrogen have negligible overlap with the plutonium 5f electrons [20], indicating that these effects could be less significant for the H/PuO 2 systems of interest here.
These issues make it difficult to determine basic kinetic properties of the initial hydriding reactions, such as hydrogen diffusion activation barriers and hopping rates. In particular, quantum nuclear vibrational effects (i.e., zero-point energy effects and/or tunneling) could prove large for these systems. Quantum zero-point energies can cause computed hydrogen surface adsorption energies in α -uranium to shift by as much as 30% [21]. Quantum tunneling has been shown to enhance hydrogen diffusion on transition metal surfaces [22] and on polycrystalline water ices by several orders of magnitude at low temperatures [23]. Recent path integral simulations have shown that quantum nuclear vibrational effects can have similar results for proton diffusion in high-pressure transition metal clathrates [24]. DFT calculations on plutonium containing materials in general are far too computationally intensive to allow for these types of calculations, which can require millions of force evaluations per diffusion pathway. Thus, plutonium hydriding studies could greatly benefit from determination of the appropriate level of DFT theory for computing hydrogen/substrate interactions, as well as a method for placing a bound on the importance of quantum nuclear vibrational effects.
In this work, we attempt to reconcile these issues by first assessing different levels of theory in order to determine whether more computationally efficient DFT calculations could yield satisfactory results for hydrogen absorbed in PuO 2 . We report on results for bulk PuO 2 from several different functionals and spin states, using GGA + U with either spin polarization or spin–orbit coupling in order to assess the effects of each. We then compute hydrogen absorption energies in the dilute limit and determine the relative importance of GGA + U and treatment of spin on hydrogen-PuO 2 interactions. Finally, we use a newly developed semi-empirical quantum Density Functional Tight Binding (DFTB) model to seed DFT calculations of the hydrogen diffusion activation energies. We use these values in a simple quantum double well formalism in order to place a bound on the quantum tunneling effects on single hydrogen diffusion in a PuO 2 lattice. Our results serve as an initial estimate for the importance of considering hydrogen quantum tunneling in these systems, which could play a significant role that requires further consideration.

2. Computational Details

PuO 2 has been determined experimentally to have a cubic lattice with F m 3 ¯ m symmetry (Figure 1). The high symmetry of the crystal yields two possible absorption sites for hydrogen: (1) a lower energy “interstitial” site, where an O–H bond is formed and is oriented towards the center of a nearby Pu-Pu axis (labeled int), and (2) a higher energy “octahedral” site, where the hydrogen ion sits between two Pu atoms (labeled oct). We focus our efforts specifically on the energetics of these interactions, as well as the the kinetics of proton hops between these sites.
All DFT calculations reported in this work were performed with the Vienna ab initio Simulation Package (VASP) [25,26,27] using projector-augmented wave function (PAW) pseudopotentials [28,29]. We have surveyed three different exchange-correlation functionals: the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) functional [30], the reformulation of PBE for solids and surfaces (PBESol) [31], and the Strongly Constrained and Appropriately Normed meta-GGA semilocal density functional (SCAN) [32]. Similar to previous studies [20], we have used the D3 dispersion energy correction [33] in all calculations discussed here, regardless of treatment of spin or choice of functional.
Lattice optimizations were performed on PuO 2 alone, using a four formula unit (12 atom) cubic unit cell with a planewave cutoff of 500 eV, Mohkhorst-Pack [34] k-point mesh of 6 × 6 × 6 , energy convergence criteria of 10 6 eV and force convergence criteria of 10 2 eV/Å, which yielded system pressures far below 1 GPa. Tests with PBE using a planewave cutoff of 600 eV, k-point mesh of 10 × 10 × 10 , energy criteria of 10 8 eV, and force criteria of 10 3 eV/Å indicated convergence of energies within 3 × 10 3 eV per formula unit and lattice constants within 0.01 Å. All lattice optimizations were performed anisotropically, i.e., without any constraints on the final symmetry of the supercell. In general, this resulted in off-diagonal lattice vector components on the order of 10 6 Å. Hence, all optimized supercells were assumed to exhibit at least orthorhombic symmetry. Hydrogen absorption calculations were performed for a 32 formula (96 atom) unit simulation cell, using a k-point mesh of 2 × 2 × 2 , energy criteria of 10 5 10 6 eV, and force criteria of 5 × 10 2 eV/Å.
Use of a linear response approach [35] to compute U eff = U J with PBE and PBESol yielded a value of 3.99 eV for both functionals. However, this value resulted in lattice constants and band gaps below experiment for all magnetic states discussed here. Hence, we have used a value of U eff = 6.00 eV for our calculations with U = 6.00 eV and J = 0.00 eV, based on results from previous work [13] and independent of choice of functional.
Given the number of different DFT calculations discussed in this work, we adopt a nomenclature to differentiate between different features. First, we consider different collinear and noncollinear periodic arrangements of spin moments, i.e., spin geometries (Table 1). Diamagnetic states (no magnetic moment) are labeled DM, ferromagnetic states are FM, and antiferromagnetic states are AF. FM states are further labeled (001) for collinear and (111) for three-dimensional noncollinear. For the AF states, we consider one-dimensional collinear geometries which we label ‘1k’, and three-dimensional noncollinear geometries that we label ‘3k’. Each of these can be further differentiated with either longitudinal (long) or transverse (trans) domains, depending on the symmetry of the spin geometry [13]. For each DFT result, we first state the type of magnetism present, DM, FM, or AF, followed by the spin configuration in parentheses, i.e., (001) or (111) for FM systems, and (1klong), (1ktrans), or (3klong) for AF systems. Thus, the following magnetic states are considered in this work: DM, FM (001), FM (111), AF (1klong), AF (1ktrans), AF (3ktrans).
We have performed calculations using both spin polarization (SP) and spin–orbit coupling (SOC). We note that in the absence of SOC the energy of the system in principle will not depend on the direction of the magnetic moment since the spin moments are not coupled to the orbital angular momenta. Hence, we would expect the total energy to be independent of choice of spin quantization axis. However, the final energy, geometry, and lattice vectors achieved during optimization can still be sensitive to choice of initial spin geometry [19], where we have observed changes in the magnitude and direction of the spin moments to occur regularly during our calculations. This effect can be pronounced when asymmetries are introduced into the lattice, such as a hydrogen point defect. As a result, we have conducted noncollinear SP calculations as a systematic way of providing additional sampling of SP spin states. The majority of our magnetic calculations have been performed at the SP level. Hence, results from calculations with spin interactions should be assumed to be from SP unless labeled with SOC. Similarly, calculations performed with GGA + U are labeled as such, and the absence of this label implies no GGA + U in that result.
DFTB calculations with self-consistent charges (SCC) were performed using the DFTB+ code [36], with a many-body repulsive energy ( E Rep ) determined from the Chebyshev Interaction Model for Efficient Simulation (ChIMES) [37,38,39,40]. Briefly, the DFTB total energy is derived from an expansion of the Kohn-Sham DFT energy to second order in charge flucuations [41]. This yields the following total energy expression:
E DFTB = E BS + E Coul + E Rep ChIMES + E Disp .
Here, E BS corresponds to the band structure energy computed from diagonalization of the DFTB Hamiltonian, and E Coul is the self-consistent charge transfer term. E Rep ChIMES is the empirical “repulsion” term determined from ChIMES, and E Disp is the dispersion energy, determined by the D3 method in our study. Additional details regarding our DFTB parameterization are discussed in the Supplementary Information.
ChIMES repulsive energies are determined through fitting linear combinations of Chebyshev polynomials to DFT reference data [38,39,40,42]. We have created an E Rep by fitting a ChIMES potential of up to two-body polynomial order of 12 and three-body polynomial order of 8 order to 75 DFT-MD snapshots for bulk PuO 2 . These were run at 600 K and at both the experimental density and under 5% compression using DM calculations and no GGA + U interaction in order to create a sufficient data set for our potential. We also have included two PuH 2 DFT-MD snapshots in order to add some sampling of Pu-H interactions, computed at 1000 K at the experimental density and under 5% compression, and also without spin or GGA + U. These minimal DFT calculations allowed us to sample a large enough part of the phase space with sufficient accuracy for the purposes of this work. The H-H, O-H, and O-O repulsive interaction were not sampled in our training set and were consequently taken from the DFTB+ miomod-hh-0-1 and mio-1-1 parameter sets. All DFTB/ChIMES calculations were subsequently performed with an SP level of theory, using periodic boundary conditions and the same k-point sampling as used for our DFT calculations. The main goal of our DFTB/ChIMES calculations was to determine initial values for diffusion hop activation energies that could then be refined by more computationally expensive DFT calculations. Further details of the ChIMES formalism and our DFTB/ChIMES model are included in the Supplementary Materials. Refinement of our DFTB/ChIMES model is the subject of future work.

3. Results

We have computed a survey of results from PBE, PBESol, and SCAN by determining optimized ground-state PuO 2 properties for a variety of spin geometries and with GGA + U (Table 2). All of our results for all functionals and levels of theory yield lattice constants that are close to the experimental value and lattice volumes that agree within ∼3%. We find that all magnetic states (FM and AF) are essentially isoenergetic for all three functionals, while the diamagnetic state exists over 3 eV per formula unit higher in energy for PBE and PBESol, and over 5.5 eV for SCAN. In addition, the FM (111) state is 1.4 × 10 2 eV higher than all AF states for PBE and PBESol, but is 8.8 × 10 3 eV lower in energy for SCAN. However, these energy differences are negligible at 298 K (~0.025 eV).
As expected, all three AF calculations with SP yield virtually isoenergetic structures, though there are differences in the lattice constants and in some cases the final lattice volume. In all cases, the AF (1klong) and AF (1ktrans) spin states resulted in a tetragonal lattice distortion of similar magnitude (0.06 Å). In contrast, the AF (3klong) state yields a smaller tetragonal distortion for PBESol (0.02 Å), an orthorhombic unit cell for PBE, and a cubic unit cell for SCAN. In general, PBESol appears to yield the best overall agreement with the experimental lattice constant and volume for the magnetic states studied here, though the differences between PBE and PBESol do not appear to be substantial. We note that U eff can be adjusted in order to yield results close to the experimental lattice volume.
Some additional differences are observed with some of the band gap results. All three AF band gaps from PBESol agree with a value of 2.92 eV. However, both PBE and SCAN yield a slightly higher band gap for AF (3klong) state relative to AF (1klong) and AF (1ktrans), with values of 2.85 vs. 2.80 eV for PBE and 2.75 vs. 2.68 eV for SCAN. PBESol yields a higher gap for the DM state relative to the AF states (3.07 eV), whereas PBE yields a lower gap (2.60 eV), and SCAN yields a metallic system. All three functionals yield a lower gap for FM (111) compared to the AF states, with a value of 2.54 eV from PBESol, 2.60 eV from PBE, and 2.16 eV for SCAN. Overall, the AF results from PBE are generally closest to the experimental value of 2.82 eV, though all three functionals yield reasonable results. Results for AF (1klong) GGA + U calculations from DFTB/ChIMES yielded an accurate lattice constant for a value of U eff = 4.08 eV, albeit for a metallic state. We found that increasing the value of U eff to as high as 12 eV did not create an insulating state. This is possibly due to either the minimal basis set used in DFTB calculations or the omission of greater than two-center terms in the DFTB Hamiltonian. Use of different confining radii and/or basis set functions could improve upon this result, though those effects are beyond our current scope of work.
Finally, in order to test the validity of our SP results, we have performed noncollinear calculations with SOC and the same value of U eff = 6.00 eV. We find that this set of calculations yields the AF (3klong) SOC state as the energetic minimum, with all three magnetic states yielding a cubic structure and lattice constant very close to that of experiment. However, once again, the AF (3ktrans) SOC and FM (111) SOC states are higher in energy by <0.01 eV, indicating that these spin geometries are essentially isoenergetic. In addition, we compute a FM (111) SOC band gap of 2.50 eV, which is very similar in value to the SP result. However, we compute AF (3klong) SOC and AF (3ktrans) SOC band gaps that are higher than the SP result by 0.17 eV (3.09 vs. 2.92 eV). The SOC relative energy and band gap results reported here differ somewhat from previous PBESol SOC results using the same value for U eff [13]. These are possibly due to a combination of the larger k-point mesh used in our study as well as the more stringent convergence criteria used in Ref. [13].
Based on these results, we choose to proceed with the PBESol functional for the rest of our calculations. PBESol provides accurate equilibrium lattice constants, particularly for the nearly cubic AF (3klong) magnetic state, as well as a band gap only slightly higher than that from experiment. Hence, it appears to provide the best combination of accuracy and computational efficiency, though we note than any of three functionals studied here would appear to provide reasonable results.
Next, we have computed hydrogen point defect formation energies for the int and oct sites using different magnetic states and levels of theory. These were calculated according to the equation:
E def = E PuO 2 / H def n E PuO 2 1 2 E H 2 ,
where the subscript ‘def’ corresponds to either the int or oct site, E PuO 2 / H def is the total int or oct defect system energy, n is the number of formula units in the PuO 2 /H system, E PuO 2 is the bulk energy per formula unit, and E H 2 is the molecular hydrogen formation energy. The goal here is to assess the effects of these parameters on interactions between the lattice and the hydrogen ion itself. For most spin states, we have fixed our simulation cell at the experimental lattice parameters, which are close the optimized values reported in Table 2. The only exceptions is the AF (1klong) calculation where the supercell lattice was held fixed at the dimensions reported in the table in order to provide direct comparison to previous H/PuO 2 collinear SP GGA + U results [16,20]. One additional AF (1klong) calculation was performed with the experimental lattice vectors with SP and without GGA + U, which is labeled as AF (1klong-iso).
Our results indicate a relative consistency in the hydrogen interstitial site formation energy, E int (Table 3). For the SP calculations, results without GGA + U yield a fairly narrow spread of E int , with a value of 0.63 eV for FM (001) at the low end and 0.71 eV for AF (3klong) at the high end. The effect of using GGA + U appears to be fairly small, with a result of 0.64 eV from AF (3klong) GGA + U. In this case, additional use of SOC has a similar magnitude of effect in the opposite direction, with a result of 0.70 eV for AF (3klong) GGA + U SOC. The tetragonal cell distortion in the AF (1klong) calculation does not appear to have a significant effect on E int , though we note that the breaking of the cubic symmetry likely results in the creation of additional interstitial sites and energies. Our DFTB/ChIMES result agrees relatively closely with the spread of results from PBESol with a value of 0.75 eV. Our results at all levels of theory for all magnetic states agree with the previously published SP GGA + U result from PBE of 0.54 eV from Ref. [16]. The exception presented here is the diamagnetic state, which yields a relatively high E int value of 0.97 eV. Our results also show some agreement with the range of PBE results of 0.48-1.19 eV from Ref. [20], though there is some ambiguity in that publication as some calculations might have been performed on the PuO 2 unit cell, where system size effects are likely significant.
Given the relative agreement between levels of theory, we choose to continue our study with the more efficient SP calculations, only. Results for the octahedral site formation energy, E oct , with results from collinear FM (001) and noncollinear AF (3klong) showing agreement within 0.08 eV, with values of 1.25 and 1.33 eV, respectively. Inclusion of GGA + U in the AF (3klong) GGA + U calculation yields a modest increase in E oct of 0.16 eV to a value of 1.49 eV. We note that our AF (3klong) GGA + U result is smaller than the collinear SP GGA + U results from PBE of 2.24 eV [16] and 2.12 eV [20]. These previous PBE results are in relatively good agreement with the DFTB/ChIMES result of 2.45 eV.
Regardless, our results indicate overall that both SOC and GGA + U likely play a relatively small role in hydrogen/PuO 2 lattice interactions. Hence, for the remainder of this work, we choose to perform calculations with SP, only and omit these two levels of theory for the sake of computational efficiency. This allows us to more efficiently explore the kinetics of hydrogen diffusion pathways as a function of magnetic state.
We now investigate the kinetics of three fundamental pathways in the H/PuO 2 system (Figure 2), namely:
  • Diffusion from one int site to its nearest neighbor int site (where an O–H bond is broken), which we label int-int.
  • Diffusion between int sites containing an O–H bond rotation (where the bond remains intact), which we label int-int2.
  • Diffusion hops from the oct site to the closest int site, which we label oct-int.
We have not computed int-oct diffusion (from an interstitial to an octahedral site) due to the large increase in energy for octahedral absorption, indicating that this pathway is likely only at extreme temperatures. In addition, oct-oct diffusion is highly unlikely to occur, given the closer proximity of low energy int sites.
We now approximate the activation energies, E A , for these pathways using a two-step approach. First, we have computed minimum energy pathways (MEP) from our DFTB/ChIMES model using the climbing image nudged elastic band (NEB) method [43]. All calculations were run within the Atomic Simulation Environment python interface [44,45], with the image dependent pair potential method [46] used to generate an initial guess for each pathway. Here, our SCC convergence was held to 2.72 × 10 4 eV ( 10 5 au) and the NEB forces were converged to 0.05 eV/Å. In addition, the simulation cell was held fixed at the experimental PuO 2 density with cubic symmetry. Decreasing the SCC convergence to 2.72 × 10 5 eV ( 10 6 au) and the NEB force convergence to 0.01 eV/Å resulted in an increase in the activation energy of 0.01 eV for the int-int2 rotation. Next, we performed constrained geometry DFT optimizations starting from the DFTB/ChIMES computed transition state geometries using PBESol with different spin states, where the hydrogen ion was held fixed and the surrounding PuO 2 ionic positions were allowed to relax. In doing so, we were able to approximate E A for several different spin states relatively efficiently at a DFT level. All PBESol optimizations were kept at the experimental PuO 2 cubic lattice, excluding AF (1klong) where the ionic coordinates and lattice vectors were scaled to the optimized tetrahedral distorted cell in order to facilitate comparison with previous DFT NEB results [20]. More exact calculations of MEPs in H/PuO 2 systems is the subject of future work.
Our results indicate a high degree of consistency between results for different spin states (Table 4). All magnetic states from PBESol conducted on a cubic simulation cell yield an E A of ~1.1 eV for the int-int transition, while the diamagnetic calculation yields a value of only 0.36 eV. DFTB/ChIMES yields a value in agreement with the PBESol magnetic results, with a value of 1.03 eV. Additional constrained optimization calculations were performed with PBE in order to examine the effect of choice of functional on our results. We find that results from both the diamagnetic and AF (3klong) states yield close agreement with results from PBESol, with values of 0.32 eV and 1.10 eV, respectively. Our result for the tetragonally deformed AF (1klong) geometry computed with PBE yields a value of 0.15 eV along the short lattice length, very close to the previously reported PBE activation energy of 0.13 eV for this transition [20].
Given the relative agreement between results from different magnetic states with a cubic simulation cell, we continue with the AF (3klong) spin geometry for PBESol, only. We find that for the int-int2 transition DFTB/ChIMES determines a relatively large E A relative to PBESol, with values of 0.36 eV vs. 0.06 eV. In contrast, we observe improved agreement for oct-int, where DFTB/ChIMES yields an activation energy of 0.27 eV compared to the value of 0.36 eV from PBESol. It is possible that the agreement of DFTB/ChIMES with PBESol could be improved with additional and more accurately computed training data, which is the subject of future work. We choose to use the PBESol AF (3klong) estimates for E A for the remainder of this work due to their likely higher accuracy than those from DFTB/ChIMES.
We estimate quantum effects on hydrogen diffusion by treating the diffusion pathways as a one-dimensional quantum double-well with a parabolic barrier, using an approach adapted from Ref. [47]. We start with the following piece-wise quadratic function:
V ( x ) = k x x 1 2 2 x < x L ( reactant region ) k b x 2 2 + V b x L x x R ( barrier region ) k x x 2 2 2 + V 0 x > x R ( product region )
Here, k is the spring constant for each well, k b is the spring constant for the barrier, x 1 is the position of the reactant minimum, x 2 is the position of the product minimum, V b is the parabolic barrier height, V 0 is the offset of the product well, and x L and x R define the reactant (left) and product (right) regions. For our work, we have chosen k = k b in order to ensure a smooth potential throughout, where k = 4 V b x 1 2 . The parameter V b is set to the approximate E A value from the AF (3klong) calculations. The potential is centered at zero for all three transitions studied here. For the symmetric int-int and int-int2 transitions, V 0 = 0 , and x 1 = 1 2 Δ x = x 2 , where Δ x is the total reaction path distance as determined from our DFTB/ChIMES NEB calculations. Furthermore, we compute the reaction/product positions as x L = 1 2 x 1 and x R = 1 2 x 2 . For the asymmetric oct-int transition, the same equations apply with the exception of V 0 = E oct E int and x 2 = 4 ( V b V 0 ) / k .
We can estimate the quantum tunneling rates by first assuming a quantum harmonic oscillator (Hermite polynomial) basis set and solving for the energy eigenstates. Care was taken to ensure convergence for each transition in terms of basis set size as well as the numerical limits of integration. Results for the int-int transition indicate small splittings between energy levels (i.e., <1 meV) (Figure 3), typical for a symmetric double well with a relatively large barrier where the particle tends to remain in a single well and tunneling transmission coefficients are small. The probability of hydrogen being found between wells becomes more pronounced as the eigenstates approach the top of the potential barrier, with the state just below the barrier yielding a splitting of slightly less than 0.01 eV. In addition, we compute a zero-point energy correction in each potential well of ~0.06 eV, and that the energy levels below the parabolic barrier tend to be spaced by ~0.13 eV.
The tunneling transmission coefficient for each eigenstate can be computed using the analytical expression for a parabolic barrier [48]:
T p a r Q M ( E ) = 1 + exp 2 π E ω b 1 , ω b = k b m .
The thermal reaction rate for hydrogen diffusion for a given pathway is then computed from the following expression:
k Q M ( T ) = 1 2 π Σ n ρ n T p a r Q M ( E n ) e β E n Σ n ρ n e β E n .
The summation is over the n computed eigenvalues E n and eigenstates ψ n of the system, and ρ n is the probability of finding the particle on the reactant side, namely:
ρ n = x r ( n ) ψ n ( x ) 2 d x .
Similar to Ref [47], we choose to define the reactant region limit for a given eigenstate, x r ( n ) , as its classical turning point. We note that the computed thermal reaction rate will have some dependence on definition of x r ( n ) , though the classical turning point is physically intuitive and consistent with well-known semi-classical approaches such as the Wentzel-Kramers-Brillouin (WKB) approximation.
Additional accuracy for our model double-well potentials could involve determination of flux-flux autocorrelation functions [49] and/or time-dependent wave-packet propagation to determine the scattering matrix [47,50]. Full determination of quantum nuclear vibrational effects in the condensed phase would require simulations using path integrals (e.g., Ref. [51]), which would necessitate the development of an accurate interatomic potential due to the extreme computational expense of these simulations [52]. Regardless, results from Ref. [47] indicate that the eigenstate transmission approach we have adopted here likely yields a high degree of accuracy for proton double-well transmission barriers as low as ~0.3 eV with respect to scattering matrix calculations.
We now make comparison to classical rate constants estimated from our approximate E A values from the PBESol AF (3klong) calculations. Classical diffusion rates were estimated from a modified version of the Eyring equation:
k Cl = k B T h exp E A / k B T .
Here, we have assume that the free energy of activation is equal to our estimate for E A and that the Eyring transmission coefficient is equal to one. This approach yields a kinetic prefactor of 6.3 × 10 12 s 1 at 300 K, which is of the same magnitude as DFT calculations for hydrogen ion diffusion on transition metal oxide surfaces (e.g., Ref. [53]).
Our results for the diffusion rates (Figure 4) indicate that int-int diffusion is likely low even at 300 K, where we compute a k QM value of 3.9 × 10 5 s 1 and a k Cl value of 2.1 × 10 6 s 1 . We predict that the rate increases by nine orders of magnitude at 600 K, with values of k QM = 2.6 × 10 4 s 1 and k Cl = 7.2 × 10 3 s 1 . These increase by an additional four orders of magnitude at 1000 K, with a values of k QM = 1.6 × 10 8 s 1 and k Cl = 6.0 × 10 7 s 1 .
In contrast, the oct-int transition is relatively frequent at 300 K due to its relatively low barrier, with values of k QM = 9.9 × 10 8 s 1 and k Cl = 5.6 × 10 6 s 1 . These values increase more modestly with increasing temperature, with values of k QM = 8.1 × 10 10 s 1 and k Cl = 1.2 × 10 10 s 1 at 600 K and at k QM = 1.0 × 10 12 s 1 and k Cl = 3.2 × 10 11 s 1 1000 K. Nonetheless, hydrogen absorption into an octrahedral defect site is unlikely at the conditions studied here, given its high formation energy of 1.5 eV.
Finally, we compute relatively rapid int-int2 O–H bond rotations for all temperatures studied here, with values of k QM = 1.1 × 10 12 s 1 and k Cl = 6.1 × 10 11 s 1 at 300 K. The temperature effect on the rate is expected to be small, given the relatively low barrier of 0.06 eV (~700 K). In this case, we compute rates of k QM = 5.3 × 10 12 s 1 and k Cl = 3.9 × 10 12 s 1 at 600 K, and k QM = 1.2 × 10 13 s 1 and k Cl = 1.0 × 10 13 s 1 at 1000 K. This transition is likely to only play a secondary role in hydrogen mass transport through the PuO 2 lattice, given that the O–H bond stays intact throughout this pathway.
The ratio of the quantum to classical rate indicate significant quantum effects for the int-int and oct-int transitions (Figure 5), both of which involve the breaking and/or forming of a covalent O–H bond. The oct-int diffusion pathways appears to show the strongest enhancement at low temperature due to quantum tunneling with a k QM / k Cl ratio of ~180, compared to the int-int ratio of ~19. However, our approach uses an identical prefactor for all pathways at all temperatures, whereas harmonic transition state theory could indicate that the prefactor for the int-int path is considerably higher than that of oct-int due to the presence of a high frequency O–H bond in the int absorption site. Regardless, our results provide reasonable estimates for the persistence of quantum nuclear vibrational effects even at high temperature. Both int-int and oct-int pathways show a quantum enhancement of approximately a factor of three at 1000 K. The relatively low E A of the int-int2 transition yields more modest quantum effects at 300 K of slightly less than 2 × enhancement, which approaches a value close to one at 1000 K. Overall, our results indicate that quantum nuclear vibrational tunneling effects could enhance hydrogen diffusion in plutonium oxide materials by more than one order of magnitude at ambient conditions. Hence, it is likely that these effects need to be included in order to make accurate predictions for experiments. More accurate quantification of these effects can accessed through the development of reactive molecular dynamics models, which is the subject of future effort.

4. Conclusions

We have used several different DFT functionals and levels of theory in order to ascertain the necessary components for accurate calculation of hydrogen diffusion in bulk PuO 2 . Our results indicate that bulk properties such as the lattice constant and band gap can have some dependence on choice of functional, level of theory (spin polarized vs. spin–orbit coupling), and magnetic state for a given value of U eff in GGA + U calculations. However, these options have much less bearing on results for H/PuO 2 interactions for systems fixed at the experimental PuO 2 volume and symmetry, where we saw only small variation in hydrogen interstitial and octahedral formation energies. This result indicates that determination of hydrogen absorption energies and diffusion parameters do not necessarily require specific spin states or more computationally intensive calculations such as those including spin–orbit coupling and/or GGA + U, which can allow for higher computational throughput approaches in the future for similar systems.
Finally, we have used approximate activation energies computed from DFT to place a quantitative bound on the effect hydrogen quantum vibrational tunneling on diffusion, which is largely ignored in such systems. Our one-dimensional models for the diffusion hop potential energy surfaces indicate that diffusion rates in H/PuO 2 systems could be enhanced by more than one order of magnitude at ambient conditions. In addition, quantum effects could prove significant even at temperatures as high as 1000 K. The method we propose here can be used as a fast screening tool for assessing possible quantum nuclear vibrational effects in any number of systems, including different hydrogen/metal or metal oxide systems and molecular crystals of low-Z organic materials. Future efforts will involve development of molecular dynamics models for H-Pu-O systems that will allow us to include anharmonic interactions, different Pu-O stoichiometries, and additional crystal defects. Our efforts will allow us to make more direct contact with experiments, where there is historically a significant reliance on quantum calculations for an atomic-level interpretation of macro-scale measurements.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/app122111005/s1, References [54,55,56,57] are cited in the supplementary materials.

Author Contributions

Conceptualization, N.G. and J.L.B.; Methodology, N.G., L.Z.-R. and R.G.M.; Software, N.G., R.K.L., C.H.P. and L.E.F.; Writing—original draft, N.G.; Writing—review & editing, N.G., L.Z.-R., R.G.M., R.K.L., C.H.P., L.E.F. and J.L.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

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. Computations were performed at LLNL using the quartz massively parallel computer. Released under document number LLNL-JRNL-840447.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
SPspin polarization
SOCspin–orbit coupling
DMdiamagentic system
FMferromagnetic system
AFantiferromagnetic system
(1klong)collinear antiferromagentic system with longitudinal geometry
(1ktrans)collinear antiferromagentic system with transverse geometry
(3ktrans)noncollinear antiferromagentic system with longitudinal geometry
int-intHydrogen diffusion between adjacent interstitial sites involving breaking of an O–H bond.
int-int2Hydrogen diffusion between adjacent interstitial sites an O–H bond rotation.
oct-intHydrogen diffusion from an octahedral to a neighboring interstitial site.

References

  1. Haschke, J.M.; Allen, T.H.; Morales, L.A. Reaction of Plutonium Dioxide with Water: Formation and Properties of PuO2+x. Science 2000, 287, 285–287. [Google Scholar] [CrossRef]
  2. Peiman, W.; Pioro, I.; Gabriel, K.; Hosseiny, M. Thermal aspects of conventional and alternative fuels. In Handbook of Generation IV Nuclear Reactors; Pioro, I.L., Ed.; Woodhead Publishing Series in Energy; Woodhead Publishing: Sawston, UK, 2016; Chapter 18; pp. 583–635. [Google Scholar]
  3. Schmidt, G.R.; Sutliff, T.J.; Dudzinski, L.A. Radioisotope Power: A Key Technology for Deep Space Exploration. In Radioisotopes; Singh, N., Ed.; IntechOpen: Rijeka, Croatia, 2011; Chapter 20. [Google Scholar]
  4. Dinh, L.; Haschke, J.; Saw, C.; Allen, P.; McLean, W. Pu2O3 and the plutonium hydriding process. J. Nucl. Mater. 2011, 408, 171–175. [Google Scholar] [CrossRef] [Green Version]
  5. Saw, C.; Haschke, J.; Allen, P.; Mclean, W.; Dinh, L. Hydrogen corrosion of plutonium: Evidence for fast grain-boundary reaction and slower intragrain reaction. J. Nucl. Mater. 2012, 429, 128–135. [Google Scholar] [CrossRef] [Green Version]
  6. Tang, J.; Qiu, R.; Chen, J.; Ao, B. Revealing the microscopic mechanism of PuO2 and α-Pu2O3 in resisting plutonium hydrogenation via ab initio molecular dynamics simulation. J. Alloy. Compd. 2021, 874, 159904. [Google Scholar] [CrossRef]
  7. Yang, Y.; Zhang, P. Hydriding and dehydriding energies of PuHx from ab initio calculations. Phys. Lett. A 2015, 379, 1649–1653. [Google Scholar] [CrossRef]
  8. Li, S.; Guo, Y.; Ye, X.; Gao, T.; Ao, B. Structural, magnetic, and dynamic properties of PuH2+x (x=0, 0.25, 0.5, 0.75, 1): A hybrid density functional study. Int. J. Hydrog. Energy 2017, 42, 30727–30737. [Google Scholar] [CrossRef]
  9. Huda, M.N.; Ray, A.K. A density functional study of atomic hydrogen adsorption on plutonium layers. Phys. B Condens. Matter 2004, 352, 5–17. [Google Scholar] [CrossRef] [Green Version]
  10. Huda, M.N.; Ray, A.K. An ab initio study of H2 interaction with the Pu (100) surface. Phys. B Condens. Matter 2005, 366, 95–109. [Google Scholar] [CrossRef]
  11. 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. [Google Scholar] [CrossRef] [Green Version]
  12. Nakamura, H.; Machida, M. Hybrid Density Functional Study on Plutonium Dioxide. In Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2013), Tokyo, Japan, 5–9 August 2013. [Google Scholar]
  13. Pegg, J.T.; Shields, A.E.; Storr, M.T.; Wills, A.S.; Scanlon, D.O.; de Leeuw, N.H. Hidden magnetic order in plutonium dioxide nuclear fuel. Phys. Chem. Chem. Phys. 2018, 20, 20943–20951. [Google Scholar] [CrossRef]
  14. 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 1998, 57, 1505. [Google Scholar] [CrossRef]
  15. Nakamura, H.; Machida, M.; Kato, M. Effects of spin–orbit coupling and strong correlation on the paramagnetic insulating state in plutonium dioxides. Phys. Rev. B 2010, 82, 155131. [Google Scholar] [CrossRef]
  16. Hernandez, S.; Holby, E.F. DFT+U Study of Chemical Impurities in PuO2. J. Phys. Chem. C 2016, 120, 13095–13102. [Google Scholar] [CrossRef]
  17. Pegg, J.T.; Shields, A.E.; Storr, M.T.; Wills, A.S.; Scanlon, D.O.; de Leeuw, N.H. Magnetic structure of UO2 and NpO2 by first-principle methods. Phys. Chem. Chem. Phys. 2019, 21, 760–771. [Google Scholar] [CrossRef] [Green Version]
  18. Dorado, B.; Amadon, B.; Freyss, M.; Bertolus, M. DFT + U calculations of the ground state and metastable states of uranium dioxide. Phys. Rev. B 2009, 79, 235125. [Google Scholar] [CrossRef]
  19. Dorado, B.; Jomard, G.; Freyss, M.; Bertolus, M. Stability of oxygen point defects in UO2 by first-principles DFT +U calculations: Occupation matrix control and Jahn-Teller distortion. Phys. Rev. B 2010, 82, 035114. [Google Scholar] [CrossRef]
  20. Zhang, L.; Sun, B.; Zhang, Q.; Liu, H.; Liu, K.; Song, H. First-Principles Study of the Hydrogen Resistance Mechanism of PuO2. ACS Omega 2020, 5, 7211–7218. [Google Scholar] [CrossRef]
  21. Soshnikov, A.; Kulkarni, A.; Goldman, N. Elucidating the Initial Steps in α-Uranium Hydriding Using First-Principles Calculations. Langmuir 2022, 38, 9335–9346. [Google Scholar] [CrossRef]
  22. McIntosh, E.M.; Wikfeldt, K.T.; Ellis, J.; Michaelides, A.; Allison, W. Quantum Effects in the Diffusion of Hydrogen on Ru(0001). J. Phys. Chem. Lett. 2013, 4, 1565–1569. [Google Scholar] [CrossRef] [Green Version]
  23. Kuwahata, K.; Hama, T.; Kouchi, A.; Watanabe, N. Signatures of Quantum-Tunneling Diffusion of Hydrogen Atoms on Water Ice at 10 K. Phys. Rev. Lett. 2015, 115, 133201. [Google Scholar] [CrossRef]
  24. Wang, H.; Yao, Y.; Peng, F.; Liu, H.; Hemley, R.J. Quantum and Classical Proton Diffusion in Superconducting Clathrate Hydrides. Phys. Rev. Lett. 2021, 126, 117002. [Google Scholar] [CrossRef] [PubMed]
  25. Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B 1993, 47, 558–561. [Google Scholar] [CrossRef] [PubMed]
  26. Kresse, G.; Hafner, J. Ab initio molecular-dynamics simulation of the liquid-metal-amorphous-semiconductor transition in germanium. Phys. Rev. B 1994, 49, 14251–14271. [Google Scholar] [CrossRef] [PubMed]
  27. Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 1996, 54, 11169–11186. [Google Scholar] [CrossRef]
  28. Blöchl, P.E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953–17979. [Google Scholar] [CrossRef] [Green Version]
  29. Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 1999, 59, 1758–1775. [Google Scholar] [CrossRef]
  30. Perdew, J.P.; Burke, K.; Enzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865–3868. [Google Scholar] [CrossRef] [Green Version]
  31. Perdew, J.P.; Ruzsinszky, A.; Csonka, G.I.; Vydrov, O.A.; Scuseria, G.E.; Constantin, L.A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett. 2008, 100, 136406. [Google Scholar] [CrossRef] [Green Version]
  32. Sun, J.; Ruzsinszky, A.; Perdew, J.P. Strongly Constrained and Appropriately Normed Semilocal Density Functional. Phys. Rev. Lett. 2015, 115, 036402. [Google Scholar] [CrossRef] [Green Version]
  33. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. [Google Scholar] [CrossRef]
  34. Monkhorst, H.J.; Pack, J.D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188–5192. [Google Scholar] [CrossRef]
  35. Cococcioni, M.; de Gironcoli, S. Linear response approach to the calculation of the effective interaction parameters in the LDA+U method. Phys. Rev. B 2005, 71, 035105. [Google Scholar] [CrossRef] [Green Version]
  36. Hourahine, B.; Aradi, B.; Blum, V.; Bonafé, F.; Buccheri, A.; Camacho, C.; Cevallos, C.; Deshaye, M.Y.; Dumitrică, T.; Dominguez, A.; et al. DFTB+, a software package for efficient approximate density functional theory based atomistic simulations. J. Chem. Phys. 2020, 152, 124101. [Google Scholar] [CrossRef] [PubMed]
  37. Lindsey, R.K.; Fried, L.E.; Goldman, N. Application of the ChIMES Force Field to Nonreactive Molecular Systems: Water at Ambient Conditions. J. Chem. Theory Comput. 2019, 15, 436–447. [Google Scholar] [CrossRef] [PubMed]
  38. Goldman, N.; Aradi, B.; Lindsey, R.K.; Fried, L.E. Development of a Multicenter Density Functional Tight Binding Model for Plutonium Surface Hydriding. J. Chem. Theory. Comput. 2018, 14, 2652–2660. [Google Scholar] [CrossRef]
  39. Lindsey, R.K.; Bastea, S.; Goldman, N.; Fried, L.E. Investigating 3,4-bis(3-nitrofurazan-4-yl)furoxan detonation with a rapidly tuned density functional tight binding model. J. Chem. Phys. 2021, 154, 164115. [Google Scholar] [CrossRef]
  40. Goldman, N.; Kweon, K.E.; Sadigh, B.; Heo, T.W.; Lindsey, R.K.; Pham, C.H.; Fried, L.E.; Aradi, B.; Holliday, K.; Jeffries, J.R.; et al. Semi-Automated Creation of Density Functional Tight Binding Models through Leveraging Chebyshev Polynomial-Based Force Fields. J. Chem. Theory Comput. 2021, 17, 4435–4448. [Google Scholar] [CrossRef]
  41. Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys. Rev. B 1998, 58, 7260–7268. [Google Scholar] [CrossRef]
  42. Pham, C.H.; Lindsey, R.K.; Fried, L.E.; Goldman, N. High-Accuracy Semiempirical Quantum Models Based on a Minimal Training Set. J. Phys. Chem. Lett. 2022, 13, 2934–2942. [Google Scholar] [CrossRef]
  43. 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. [Google Scholar] [CrossRef]
  44. Larsen, A.H.; Mortensen, J.J.; Blomqvist, J.; Castelli, I.E.; Christensen, R.; Dulak, M.; Friis, J.; Groves, M.N.; Hammer, B.; Hargus, C.; et al. The atomic simulation environment—A Python library for working with atoms. J. Phys. Condens. Matter 2017, 29, 273002. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Bahn, S.R.; Jacobsen, K.W. An object-oriented scripting interface to a legacy electronic structure code. Comput. Sci. Eng. 2002, 4, 56–66. [Google Scholar] [CrossRef] [Green Version]
  46. Smidstrup, S.; Pederson, A.; Stokbro, K.; Jónsson, H. Improved initial guess for minimum energy path calculations. J. Chem. Phys. 2014, 140, 214106. [Google Scholar] [CrossRef] [Green Version]
  47. Garashchuk, S.; Gu, B.; Mazzuca, J. Calculation of the Quantum-Mechanical Tunneling in Bound Potentials. J. Theor. Chem. 2014, 2014, 240491. [Google Scholar] [CrossRef] [Green Version]
  48. Landau, L.D.; Lifshitz, E.M. Quantum Mechanics: Non-Relativistic Theory; Pergamon Press: Oxford, UK, 1977. [Google Scholar]
  49. Miller, W.H.; Schwartz, S.D.; Tromp, J.W. Quantum Mechanical Rate Constants for Bimolecular Reactions. J. Chem. Phys. 1983, 79, 4889–4898. [Google Scholar] [CrossRef]
  50. Trybel, F.; Cosacchi, M.; Meier, T.; Axt, V.M.; Steinle-Neumann, G. Proton dynamics in high-pressure ice-VII from density functional theory. Phys. Rev. B 2020, 102, 184310. [Google Scholar] [CrossRef]
  51. Glaesemann, K.R.; Fried, L.E. Quantitative molecular thermochemistry based on path integrals. J. Chem. Phys. 2005, 123, 034103. [Google Scholar] [CrossRef]
  52. Kimizuka, H.; Thomsen, B.; Shiga, M. Artificial neural network-based path integral simulations of hydrogen isotope diffusion in palladium. J. Phys. Energy 2022, 4, 034004. [Google Scholar] [CrossRef]
  53. Li, S.C.; Zhang, Z.; Sheppard, D.; Kay, B.D.; White, J.M.; Du, Y.; Lyubinetsky, I.; Henkelman, G.; Dohnálek, Z. Intrinsic Diffusion of Hydrogen on Rutile TiO2(110). J. Am. Chem. Soc. 2008, 130, 9080–9088. [Google Scholar] [CrossRef]
  54. Wang, Y.; Shepler, B.C.; Braams, B.J.; Bowman, J.M. Full-dimensional, ab initio potential energy and dipole moment surfaces for water. J. Chem. Phys. 2009, 131, 054511. [Google Scholar] [CrossRef]
  55. Wang, Y.; Huang, X.; Shepler, B.C.; Braams, B.J.; Bowman, J.M. Flexible, ab initio potential, and dipole moment surfaces for water. I. Tests and applications for clusters up to the 22-mer. J. Chem. Phys. 2011, 134, 094509. [Google Scholar] [CrossRef] [PubMed]
  56. Lindsey, R.K.; Fried, L.E.; Goldman, N. ChIMES: A Force Matched Potential with Explicit Three-Body Interactions for Molten Carbon. J. Chem. Theory Comput. 2017, 13, 6222–6229. [Google Scholar] [CrossRef] [PubMed]
  57. Tersoff, J. Empirical interatomic potential for carbon, with application to amorphous-carbon. Phys. Rev. Lett. 1988, 61, 2879. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Snapshot of a 32 formula unit PuO 2 supercell, with zoomed in images of the interstitial (int) and octahedral (oct) hydrogen bulk absorption sites. Periodic boundaries for bulk PuO 2 are shown by the grey lines. Pu atoms are colored grey, oxygen atoms are red, and the hydrogen atom is white.
Figure 1. Snapshot of a 32 formula unit PuO 2 supercell, with zoomed in images of the interstitial (int) and octahedral (oct) hydrogen bulk absorption sites. Periodic boundaries for bulk PuO 2 are shown by the grey lines. Pu atoms are colored grey, oxygen atoms are red, and the hydrogen atom is white.
Applsci 12 11005 g001
Figure 2. Hydrogen diffusion pathways investigated in this study. The hydrogen in the octahedral site in (c) is colored blue for the sake of clarity.
Figure 2. Hydrogen diffusion pathways investigated in this study. The hydrogen in the octahedral site in (c) is colored blue for the sake of clarity.
Applsci 12 11005 g002
Figure 3. Hydrogen probability densities for the int-int diffusion hop. Probability densities are shifted upward by their eigenvalue and colored for the sake of clarity. The double well results in a splitting of all energy levels, which is on a scale smaller scale than shown in the plot for lower eigenstates.
Figure 3. Hydrogen probability densities for the int-int diffusion hop. Probability densities are shifted upward by their eigenvalue and colored for the sake of clarity. The double well results in a splitting of all energy levels, which is on a scale smaller scale than shown in the plot for lower eigenstates.
Applsci 12 11005 g003
Figure 4. Hydrogen diffusion rates for the three pathways investigated in this study. Results for int-int are in black, int-int2 in red, and oct-int in blue. Solid lines denote the quantum result and dashed lines the classical result.
Figure 4. Hydrogen diffusion rates for the three pathways investigated in this study. Results for int-int are in black, int-int2 in red, and oct-int in blue. Solid lines denote the quantum result and dashed lines the classical result.
Applsci 12 11005 g004
Figure 5. Ratio of the hydrogen quantum to classical diffusion rates for the three pathways investigated in this study. Results for int-int and oct-int diffusion hops are shown on a log scale for the sake of clarity.
Figure 5. Ratio of the hydrogen quantum to classical diffusion rates for the three pathways investigated in this study. Results for int-int and oct-int diffusion hops are shown on a log scale for the sake of clarity.
Applsci 12 11005 g005
Table 1. Magnetic states investigated in this study. Parenthetical terms correspond to magnetic orderings defined in the text.
Table 1. Magnetic states investigated in this study. Parenthetical terms correspond to magnetic orderings defined in the text.
Ion (Fractional Coord.)FM 001 FM 111 AF (1klong)AF (3klong)AF (1ktrans)AF (3ktrans)
0 , 0 , 0 0 , 0 , 1 1 , 1 , 1 0 , 0 , 1 1 , 1 , 1 0 , 0 , 1 1 , 1 , 1
1 2 , 1 2 , 0 0 , 0 , 1 1 , 1 , 1 0 , 0 , 1 1 ¯ , 1 ¯ , 1 0 , 0 , 1 ¯ 1 , 1 ¯ , 1 ¯
1 2 , 0 , 1 2 0 , 0 , 1 1 , 1 , 1 0 , 0 , 1 ¯ 1 ¯ , 1 , 1 ¯ 0 , 0 , 1 1 ¯ , 1 ¯ , 1
0 , 1 2 , 1 2 0 , 0 , 1 1 , 1 , 1 0 , 0 , 1 ¯ 1 , 1 ¯ , 1 ¯ 0 , 0 , 1 ¯ 1 ¯ , 1 , 1 ¯
Table 2. Physical properties of bulk PuO 2 as a function of DFT functional and magnetic state. All DFT results reported here include GGA + U with U eff = 6.00 eV. Calculations were performed with SP, excluding those denoted with SOC. Listing of a single lattice constant value corresponds to optimization to a cubic lattice within an error of 0.01 Å per lattice length. The DFTB/ChIMES result was computed with U eff = 4.08 eV in order to yield an accurate lattice constant.
Table 2. Physical properties of bulk PuO 2 as a function of DFT functional and magnetic state. All DFT results reported here include GGA + U with U eff = 6.00 eV. Calculations were performed with SP, excluding those denoted with SOC. Listing of a single lattice constant value corresponds to optimization to a cubic lattice within an error of 0.01 Å per lattice length. The DFTB/ChIMES result was computed with U eff = 4.08 eV in order to yield an accurate lattice constant.
E XC Magnetic StateLattice Const (Å)Volume (Å 3 )Band Gap (eV)Relative Energy to AF (3klong) (eV/F.U.)
PBESol with GGA + UDM5.39, 5.35, 5.37154.853.073.06
FM (111)5.41, 5.39, 5.37156.592.54 1.4 × 10 2
AF (1klong)5.37, 5.37, 5.43156.582.92 7.5 × 10 4
AF (3klong)5.38, 5.40, 5.38156.302.92
AF (1ktrans)5.37, 5.43, 5.37156.582.92 7.5 × 10 4
PBESol SOC with GGA + UFM (111)5.39156.592.50 7.5 × 10 3
AF (3klong)5.39156.593.09
AF (3ktrans)5.39156.593.09 5.0 × 10 4
PBE with GGA + UDM5.43160.102.603.06
FM (111)5.45161.882.60 1.4 × 10 2
AF (1klong)5.42, 5.42, 5.48160.982.80 1.8 × 10 3
AF (3klong)5.42, 5.47, 5.44161.282.85
AF (1ktrans)5.42, 5.48, 5.42160.982.80 1.8 × 10 3
SCAN with GGA + UDM5.39156.590.05.52
FM (111)5.40157.462.16 8.8 × 10 3
AF (1klong)5.41, 5.41, 5.47160.102.68 3.4 × 10 2
AF (3klong)5.43160.102.75
AF (1ktrans)5.41, 5.47, 5.41160.102.67 3.4 × 10 2
DFTB/ChIMES with GGA + UAF (1klong)5.39156.590.0
Experiment [1]5.396157.032.80
Table 3. Hydrogen interstitial and octahedral site formation energies from PBESol. Calculations were performed with SP only and exclude GGA + U and SOC unless otherwise noted. The AF (1klong) calculations was performed on the tetargonally distorted lattice, whereas the AF (1klong-iso) calculation was performed with the cubic experimental lattice.
Table 3. Hydrogen interstitial and octahedral site formation energies from PBESol. Calculations were performed with SP only and exclude GGA + U and SOC unless otherwise noted. The AF (1klong) calculations was performed on the tetargonally distorted lattice, whereas the AF (1klong-iso) calculation was performed with the cubic experimental lattice.
Defect Magnetic StateE (eV)
InterstitialDiamagnetic0.97
FM (001)0.63
FM (111)0.71
AF (1klong)0.64
AF (1klong-iso)0.61
AF (3klong)0.71
AF (3klong) GGA + U0.64
AF (3klong) GGA + U SOC0.70
DFTB/ChIMES0.75
OctahedralFM (001)1.25
AF (3klong)1.33
AF (3klong) GGA + U1.49
DFTB/ChIMES2.45
Table 4. Estimated activation energies for hydrogen diffusion hops from spin polarization calculations, only, where both GGA + U and SOC have been omitted. All calculations were performed on the experimental cubic lattice, excluding the PBE AF (1klong) results, where the tetragonally distorted cell was used and the hydrogen was moved in the direction of one of the shorter lattice vectors.
Table 4. Estimated activation energies for hydrogen diffusion hops from spin polarization calculations, only, where both GGA + U and SOC have been omitted. All calculations were performed on the experimental cubic lattice, excluding the PBE AF (1klong) results, where the tetragonally distorted cell was used and the hydrogen was moved in the direction of one of the shorter lattice vectors.
Pathway Model Magnetic State E A (eV)
int-intDFTB/ChIMESAF (1klong)1.03
PBESolDiamagnetic0.36
FM (001)1.17
FM (111)1.10
AF (3klong)1.10
PBEDiamagentic0.32
AF (1klong)0.15
AF (3klong)1.10
int-int2DFTB/ChIMESAF (1klong)0.36
PBESolAF (3klong)0.06
oct-intDFTB/ChIMESAF (1klong)0.27
PBESolAF (3klong)0.36
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Goldman, N.; Zepeda-Ruiz, L.; Mullen, R.G.; Lindsey, R.K.; Pham, C.H.; Fried, L.E.; Belof, J.L. Estimates of Quantum Tunneling Effects for Hydrogen Diffusion in PuO2. Appl. Sci. 2022, 12, 11005. https://doi.org/10.3390/app122111005

AMA Style

Goldman N, Zepeda-Ruiz L, Mullen RG, Lindsey RK, Pham CH, Fried LE, Belof JL. Estimates of Quantum Tunneling Effects for Hydrogen Diffusion in PuO2. Applied Sciences. 2022; 12(21):11005. https://doi.org/10.3390/app122111005

Chicago/Turabian Style

Goldman, Nir, Luis Zepeda-Ruiz, Ryan G. Mullen, Rebecca K. Lindsey, C. Huy Pham, Laurence E. Fried, and Jonathan L. Belof. 2022. "Estimates of Quantum Tunneling Effects for Hydrogen Diffusion in PuO2" Applied Sciences 12, no. 21: 11005. https://doi.org/10.3390/app122111005

APA Style

Goldman, N., Zepeda-Ruiz, L., Mullen, R. G., Lindsey, R. K., Pham, C. H., Fried, L. E., & Belof, J. L. (2022). Estimates of Quantum Tunneling Effects for Hydrogen Diffusion in PuO2. Applied Sciences, 12(21), 11005. https://doi.org/10.3390/app122111005

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop