Development of a Multi-center Density Functional Tight Binding Model

Subscriber access provided by Kaohsiung Medical University

Condensed Matter, Interfaces, and Materials

Development of a Multi-center Density Functional Tight Binding Model for Plutonium Surface Hydriding Nir Goldman, Bálint Aradi, Rebecca Kristine Lindsey, and Laurence E. Fried J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00165 • Publication Date (Web): 03 Apr 2018 Downloaded from http://pubs.acs.org on April 5, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Development of a Multi-center Density Functional Tight Binding Model for Plutonium Surface Hydriding Nir Goldman,∗,†,‡ Bálint Aradi,¶ Rebecca K. Lindsey,† and Laurence E. Fried† †Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California 94550, United States ‡Department of Chemical Engineering, University of California, Davis, California 95616, United States ¶Bremen Center for Computational Materials Science, Universität Bremen, P.O.B. 330440, D-28334, Bremen, Germany E-mail: [email protected] Phone: +1-925-422-3994

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract

We detail the creation of a multi-center density functional tight binding (DFTB) model for hydrogen on δ-plutonium, using a framework of new Slater-Koster interaction parameters and a repulsive energy based on the Chebyshev Interaction Model for Efficient Simulation (ChIMES), where two and three-center atomic interactions are represented by linear combinations of Chebyshev polynomials. We find that our DFTB/ChIMES model yields a total electron density of states for bulk δ-Pu that compares well to that from Density Functional Theory, as well as to a grid of energy calculations representing approximate H2 dissociation paths on the δ-Pu (100) surface. We then perform molecular dynamics simulations and minimum energy pathway calculations to determine the energetics of surface dissociation and sub-surface diffusion on the (100) and (111) surfaces. Our approach allows for the efficient creation of multi-center repulsive energies with a relatively small investment in initial DFT calculations. Our efforts are particularly pertinent to studies that rely on quantum calculations for interpretation and validation, such as experimental determination of chemical reactivity both on surfaces and in condensed phases.

2


Page 2 of 33

Page 3 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

Introduction

Plutonium metal is known to have six allotropes between ambient and melting conditions, 1 which yields a wide range of different physical and chemical properties. The metastable facecentered-cubic (fcc) δ-phase is of particular interest for engineering applications in part due to its high ductility and malleability. 2 This phase is stabilized by mixing with small concentrations of group III metals 3 and exhibits a number of unconventional properties, including a negative coefficient of thermal expansion and the ability to become superconducting when alloyed. 4 However, δ-Pu is highly sensitive to corrosion from adsorbed atmospheric gases, which presents significant challenges due to the ubiquity of these compounds in experimental setups. Surface hydrogenation in particular is problematic, 5,6 where the Pu-hydride product can form spontaneously, 7,8 which can accelerate strongly exothermic oxidation reactions and even spontaneous ignition. 9 Despite these issues, relatively little is known about pertinent surface chemistry, and experiments would greatly benefit from detailed atomistic knowledge of the chemical reactivity of hydrogen on the surface of these systems. 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 calculations, though, are extremely computationally intensive, and are generally limited to exceedingly small system sizes and time scales. In particular, f-electron solids such as plutonium can be difficult to model efficiently, 10–17 where the high computational effort due to electron correlations and minimizing to specific electronic spin-states (e.g., anti-ferromagnetic vs. ferromagnetic) precludes the calculation of time-dependent dynamics, reaction barriers, and other relevant properties. Our previous DFT calculations for hydrogen surface reactivity on the δ-Pu (100) and (111) surfaces yielded insight into plutonium corrosion, 18 but were somewhat limited in scope due to these extreme computational costs. Purely empirical classical molecular dynamics force fields, where atomic forces are computed by parameterized potential energy surfaces, are usually highly computationally efficient but tend to yield poor results outside of their fitting regime, and do not allow for 3



the simultaneous determination of electronic properties. Thus, there exists a strong need for atomistic simulation approaches that can achieve acceptable breadth and accuracy while exhibiting a high degree of computational efficiency. To this end, we have created a density functional tight binding (DFTB) semi-empirical quantum approach for the hydrogen/δ-Pu system. DFTB utilizes a balance of empiricism with approximate quantum mechanics to yield computational times that are approximately three orders of magnitude more efficient than Kohn-Sham DFT while retaining most of its accuracy. 19,20 The formalism for DFTB with self-consistent charges (SCC) has been discussed in detail elsewhere. 19,21–25 Briefly, the method assumes neutral, spherically symmetric charge densities on the atoms and expands the DFT Hamiltonian to second-order in charge fluctuations. The DFTB total energy is expressed as EDFTB = EBS + ECoul + Erep . The first term, EBS , corresponds to the band structure energy computed via sum over occupied electronic states from the approximate DFTB Hamiltonian. EBS is usually computed from pre-tabulated Slater-Koster tables derived from calculations with a minimal basis set, where both the electronic wave functions and electron density are subjected to separate confining or compression potentials. The compression potentials force the wavefunction/electron density to zero at relatively large distances from the nuclei, which has been shown to improve transferability of the Slater-Koster tabulations. 26 The second term, ECoul , corresponds to a charge fluctuation term, and the third term, Erep , (the repulsive energy) corresponds to ion-ion repulsions, as well as Hartree and exchange-correlation double counting terms. In practice, Erep is expressed as a short-ranged empirical function whose parameters are fit to reproduce DFT or experimental data, and can be either pair-wise 25,27 or contain multi-center interactions. 28,29 In this work, we first detail our parameterization of the Slater-Koster interactions for plutonium, including the radial and density compression radii, and charge transfer and other terms to account for f-electron correlations. We discuss the determination of a multi-center empirical repulsive energy by force matching linear combinations of Chebyshev polynomi-

4


Page 4 of 33

Page 5 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


als 30,31 to our previously computed DFT results on hydrogen/δ-Pu. The use of a functional form with purely linear dependencies allows for optimal parameters to be determined in a single step, making optimizations for any number of systems extremely rapid. We then compute molecular dynamics trajectories and minimum energy pathways of molecular and dissociated hydrogen on the (100) and (111) δ-Pu surfaces. We find that our approach provides a straightforward method for including additional multi-center effects in DFTB for systems such as actinides where these interactions can be significant, as well as a general method for determining DFTB models for gas/metal corrosion for any number of systems.

2

Methods

Our previous DFT calculations calculations on the hydrogen/δ-Pu system 18 were conducted with projector-augmented wave function (PAW) pseudopotentials 32,33 and the Perdew, Burke, and Ernzerhof exchange correlation functional (PBE). 34 In short, these calculations were performed with the VASP code 35–37 utilizing a cutoff of 400 eV, fourth order Methfessel-Paxton smearing 38 with a value of 0.3 eV, and a wavefunction convergence criteria of 10−4 eV. A Monkhorst-Pack k-point mesh 39 of 5 × 5 × 5 was used for calculations on a 32 atom bulk system. The underestimation of 5f electron correlations was mitigated through use of onsite orbital Coulombic repulsions (e.g., LDA+U). 14 In this case, we used the GGA+U correction of Dudarev et al., 40 with a U − J value of 1.30 eV for the f-orbital interactions. The GGA+U approach was found necessary in order to yield a cell volume as close as possible to the experimental result. 18,41 This was necessary in order to get an accurate surface configuration which can have a strong effect on the resulting surface chemistry. Care was taken to minimize each geometry to as close to an anti-ferromagnetic surface (AFM) as possible. 12,13,42 The magnetism of δ-Pu is complicated, 43 and has been the subject of a number of recent experimental and theoretical studies. 15,16,44 Optimal parameters for the DFTB plutonium Slater-Koster interactions were chosen by

5



comparison to the total electronic density of states (DOS) from DFT computed for the optimized δ-Pu bulk system. We performed a semi-exhaustive search where we systematically varied the compression radii for the wave functions (rψα , α = s, p, d, orf ), the separate compression radii for the electronic density (rn ), and the 7f–6d charge transfer in the Pu atom (focc − docc ). All of our parameterizations used quadratic compression potentials, and relativistic effects were included using the zeroth order regular approximation (ZORA) to the Dirac equation. In this work, the s, p, and d-orbital compression was varied from 3.50 to 10.00 au, and the f-orbital compression from 3.00 to 9.50 au. The density compression for all orbitals was kept the same, and was varied from 3.00 to 9.50 au. Similar to recent parameterizations for cerium, 45 we found that changing the electronic reference state for the Pu atoms improved our predictions of the electronic density of states. Hence, we varied the d-orbital partial charge to values up to 0.5. This is consistent with a number of DFT results that have indicated non-integer orbital occupations for transition 46 and f-electron metals. 45 Lastly, we have used the pseudo self-interaction +U correction to account for the strong electron correlations in higher angular momentum orbitals, 47 where we also found it necessary to apply the same value of the +U correction to both 6d and 5f-orbitals. 45 The Pu Slater orbital basis set included the following six exponential coefficients: 1.000, 3.114, 9.695, 30.189, 94.000, 292.691. We found that use of different values (e.g., 1.000, 3.950, 15.620, 61.740, 244.000) had an insignificant effect on the resulting DOS for the bulk system. Hydrogen Slater-Koster parameters were taken from the mio-1-1 parameter set (available for download from http://dftb.org). Our DFTB self-consistent charge calculations were performed with the DFTB+ code, 23 with non-collinear spin-polarization. 48,49 We have chosen to proceed with the following orbital parameters: rψs,p,d = 4.5 au, rψf = 4.0 au, rn = 15.0 au, focc = 5.6, and docc = 0.4. The charge transfer value is consistent with a total occupied d-electron population of ∼10%, determined by integration of our DFT-computed projected electronic density of states up to the Fermi energy. We used a value of U = 0.544 eV for the pseudo self-interaction correction, which was

6


Page 6 of 33

Page 7 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tuned by hand to yield an optimal DOS result compared to DFT+U in terms of peak heights and shape of the distribution over a range of −3 to 2 eV. The Pu geometry used for these calculations was taken from our optimized AFM spin-polarized PBE+U DFT calculation so as to avoid creating a separate Erep for each set of parameters. We note that DFTB-SCC has been shown to be more sensitive to the +U correction than DFT calculations for a given system. 45 All of our DFTB calculations were minimized to the AFM surface. Comparison of the total electronic DOS for δ-Pu to that from DFT (Figure 1), indicates that our DFTB model yields a reasonably accurate representation of the electronic states. Here, we have computed the DOS with a 6 × 6 × 6 Monkhorst-Pack mesh and an SCC convergence criteria of 2.72 × 10−5 eV (10−6 au). In particular, we observe a similar peak height between the two methods at slightly less than −1 eV, and similar shoulders from −1 to approximately 2 eV. The DOS from DFTB has somewhat lower values at energies below −1 eV. However, these states are unlikely to contribute to surface chemistry. We then determined the repulsive energy for the hydrogen-plutonium system through force matching to the forces from our previous DFT calculations. The repulsive energy in this work uses the recently developed classical force field referred to as the Chebyshev Interaction Model for Efficient Simulation (“ChIMES”). 30,31 The ChIMES model represents interatomic forces by projecting the many-body DFT energy onto linear combinations of Chebyshev polynomials representing two and three-body (atomic center) interactions. Use of Chebyshev polynomials of the first kind contains a number of advantages over other approaches (e.g., a standard power series) in that: (1) one can approach basis set completeness due to their recursive nature and orthogonality, (2) their derivatives can be determined efficiently through Chebyshev polynomials of the second kind, 50 and (3) their values are bounded in the range [−1, 1] which can allow for improved conditioning of the fitting matrix. The ChIMES total energy of interaction is described as follows:

EChIMES =

XNX i

Eij +

j>i

N X XX i

7

j>i k>j


Eijk .

(1)


Page 8 of 33

Here, N is the total number of atoms in the system, Eij is the pairwise interaction energy, and Eijk is the energy between triplets of atoms. The two-body energy is expressed as follows: Eij = fpαi βj (rij ) + fcαi βj (rij )

O2 X

Cnαi βj Tn (sij )

(2)

n=1

In this case, Tn (sij ) represents a Chebyshev polynomial of order n, sij is the pair distance α βj

transformed over the interval [−1, 1] (discussed below), and Cn i

is the corresponding co-

efficient for the interaction between atom types αi and βj . Permutational invariance of the αβ

polynomials is enforced for all interactions, e.g., Cnαβ ≡ Cnβα . The term fc i j (rij ) is a smooth cutoff function which is set to zero beyond a maximum distance defined for a given atom pair αβ

max 3 type (e.g., H-Pu,), fc i j (rij ) = (1 − rij /rαβ ) . In order to prevent sampling of rij distances αβ

below what is present in our DFT training set, we introduce a penalty function fp i j (rij ), which we define as follows: 30

α βj

fp i

rpαi βj

αβ

= ps (rp i j )3

  rmin + pd − rij , rij − pd < rmin αi βj αi βj =   0, otherwise.

(3)

(4)

The parameter ps is a penalty function scaling factor, set to a value of 106 , and pd is the penalty distance, set to a value of 10−2 for this work. These value allow for atoms to be “pushed” to larger distances to avoid unphysical regions of the potential. We note that the penalty function was not sampled for any of the geometry optimization, NEB, or MD calculations presented here. We map the interatomic distances over the interval of [−1,1] by first transforming the internuclear distance rij to the Morse variable, xij = exp(−rij /λαβ ), where λ is the Morse variable range parameter, 51–53 defined individually for each type of atomic pair interaction. The Morse variables lead to a natural decrease in the interaction strength as distance is increased. For this work, we set the value of λ to 0.70 Å for H-H interactions, 3.20 Å for Pu8


Page 9 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Pu interactions, and 1.85 Å for H-Pu interactions. These correspond approximately to the nearest neighbor distance for each atom pair type from our fitting set, though we found that the results of our fit were relatively insensitive to these values. We then define the variable αβ sij to be within the range [−1, 1] through the operation sij ≡ (xij − xαβ avg )/xdiff where:

αβ αβ xαβ avg = 0.5(xmax + xmin )

(5)

αβ αβ xαβ diff = 0.5(xmax − xmin ).

(6)

min xαβ max = exp(−rαβ /λαβ )

(7)

max xαβ min = exp(−rαβ /λαβ )

(8)

Similar to the two-body representation, the three-body energy is given as the product of Chebyshev polynomials for each constituent set of atomic pairs:

Eijk =

fcαi βj

(rij ) fcαi γk

(rik ) fcβj γk

(rjk )

O3 X O3 X O3 X

0 αi βj γk Cmpq Tm (sij ) Tp (sik ) Tq (sjk ) .

(9)

m=0 p=0 q=0

The single sum given for the two-body energy is now replaced with a triple sum for the ij, ik, and jk polynomials, yielding a single permutationally invariant coefficient for each αβγ set of powers and atom types, Cmpq . The primed sum indicates that only terms for which

two or more of the m, p, q indicies are greater than zero are included, which guarantees that three bodies i, j, k enter into the expression. The expression for Eijk also contains the fc smoothly varying cutoff functions for each constituent pair distance. Penalty functions are not included for three-body interactions and are instead handled by the two-body EChIMES terms. Extension of the ChIMES formalism to include four-body interactions is the subject of future work. We determine optimal ChIMES Erep parameters for our chosen Slater-Koster interactions with a procedure similar to that used for previous DFTB force matched models. 25 First, the

9



corresponding DFTB forces from the band structure and electrostatic energies for the training set were computed using the DFTB+ code 23 with the repulsive energies set to zero. The difference from the DFT total atomic forces in each direction were then used to determine the force components for the repulsive energy, Frep . This allows for the ChIMES Erep to be determined through the overdetermined matrix equation M C = Frep . The matrix M corresponds to derivatives of the forces with respect to the coefficients, which involves the Chebyshev polynomials Tn in Eq. 2 and 9 and also their derivatives, and C corresponds to the vector of Chebyshev polynomial coefficients. The main advantage of the ChIMES formalism over other many-bond functional forms such as bond-order force fields (e.g., Ref. 54 ) is that its functional form is linear in terms of the underlying polynomial basis set. As a result, optimal parameters are solved for extremely rapidly (e.g., within minutes) through solving a singular value decomposition (SVD), as opposed to iterative non-linear least squares solvers (e.g., Levenberg-Marquardt) or global minimum searches (e.g., simulated annealing), which can take many orders of magnitude longer to find optimal sets of parameters. The SVD also allows the systematic exclusion of nearly degenerate parameter combinations through inspection of the singular values, which can improve the numerical stability of the fit. This rapidity with which a multi-center ChIMES Erep can be determined allows for an efficient method to include greater than two-body interactions in DFTB calculations, which can be poorly represented in the two-center Hamiltonian and many Erep models. For these systems, our DFT training data included “energy grid” calculations of H2 dissociation on the (100) and (111) surfaces, where the hydrogen ions were separately translated in the transverse directions and had their positions optimized relative to the surface normal, as well as to H2 molecular adsorption energies computed over a similar grid. In our training set, we used 55 dissociation configurations from the (100) surface and 135 from the (111) surface, as well as molecular adsorption calculations at a 55 different grid points on the (111) surface. This yielded a total of 235 configurations, which encompassed the bulk of our previous DFT results. The (100) surface consisted of a 32 atom slab with cell dimensions of

10


Page 10 of 33

Page 11 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


39.161 Å × 9.617 Å × 9.161 Å, with a 30 Å vacuum placed above the cell direction (Figure 2(a)). The (111) surface contained a supercell of 24 plutonium atoms, with cell dimensions of 5.895 Å × 46.232 Å × 6.379 Å, with a 30 Å vacuum similarly placed above the cell direction (Figure 2(b)). Our data set of (111) was somewhat larger due to the fact that this surface had a smaller number of atoms and was subsequently more efficient to compute. DFTB calculations for Frep were performed with the same fourth-order Methfessel-Paxton smearing of 0.3 eV, a k-point mesh of 2×6×6, and an SCC convergence criteria of 2.72×10−5 eV (10−6 au). Here, we used the mio-1-1 H-H Erep as well, which has proven to be accurate over a broad range of pressures and temperatures. 29 Hence, our ChIMES Erep was fit to H-Pu min and Pu-Pu interactions, only. The value of rαβ for each pair type was set by the minimum min min distance sampled in our data set, i.e., rHP u = 1.85, and rP uP u = 3.0. Maximum (cutoff) max pairwise distances rαβ were set to 4.5 Å for both H-Pu and Pu-Pu interactions (the second

nearest neighbor distance in the Pu lattice). Extending the radial range of our Erep did not have a noticeable effect on the accuracy of our results. Similar to our DFT calculations, a “frozen lattice” approximation was used where all Pu atoms were held fixed in their DFT optimized slab positions for all energy and optimization calculations, and only the hydrogen atoms were allowed to move. As a result, the forces acting on the hydrogen atoms were given a weight in our SVD calculation that was 104 times greater than that of the plutonium atoms. We found an optimal fit by using 16th order polynomial interactions for the two-body potential and 10th order for the three-body potential. This yielded a total of 48 two-body parameters (16 for each pair type) and 741 unique three-body parameters. This combination was found to yield the highest possible degree of accuracy while minimizing overfitting, which resulted in unstable calculations. For larger data sets, cross-validation methods can be used to determine optimal complexity in the potential energy surface. 55,56 The computational cost of our Erep calculation is insignificant due to the high cost of matrix diagonalization

11



inherent in DFTB caculations. For our ChIMES Erep , we found a root mean square (RMS) error in the forces acting on the hydrogens of 0.15 eV/Å (3.50 kcal/mol·Å). For the sake of comparison, we have computed an additional two-body only ChIMES Erep with 16th order polynomial interactions, which yielded an RMS error in the hydrogen forces of 0.35 eV/Å (7.97 kcal/mol·Å). Our “two-plus-three-body” and “two-body only” Erep models yielded RMS error for the plutonium forces of 2.05 eV/Å (47.17 kcal/mol·Å) and 1.99 eV/Å (45.82 kcal/mol·Å), respectively, though these poor values are expected because Pu-Pu interactions were largely excluded from our training set due to the frozen lattice approximation.

3

Results and Discussion

We have tested the accuracy of our two-plus-three-body ChIMES/DFTB model by first computing the energetics of H2 interactions on the δ-Pu (100) system. We find that the H2 molecule preferentially adsorbs in the center site with an energy [Eads = Eslab+adsorbate − (Eslab + Eadsorbate )] of +0.02 eV (slightly endothermic), compared to a value from DFT of −0.04 eV. Similar to DFT, the optimized H2 geometry is normal to the surface, and has a center of mass height relative to the surface of 3.18 Å, compared to the DFT value of 3.24 Å. We then computed an identical energy grid to our previous DFT calculations: we keep the bottom-most H atom positioned within the (100) center site, translate the top-most H atom along a grid of points spaced by 0.3 to 1.0 Å in each transverse direction, and then optimize both H positions relative to the surface normal, only. Our results indicate a high degree of accuracy (Figure 3), where our model yields an RMS error in the grid energies of only 0.08 eV, and an RMS error in the hydrogen positions (height above the surface) of only 0.09 Å. Calculation of the same grid of energy points using our two-body only ChIMES Erep model yielded relatively high RMS errors of 0.32 eV in the energies and 0.13 Å in the surface normal positions. Hence, we report results from our two-plus-three-body ChIMES/DFTB model only for the remainder of this work.

12


Page 12 of 33

Page 13 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We have computed the adsorption of a hydrogen ion in common adsorption sites on both (100) and (111) surfaces (Table 1). In general, our results for the adsorption energies and hydrogen ion distances from the Pu surface agree well with the spread of numbers from previous spin polarization and spin-orbit coupling DFT calculations. 5,18 We notice some variability in the total system energy based on the initial choice of AFM spin geometry. In addition, our DFTB calculations have used a stricter convergence criteria for the forces and energies as well as a larger k-point sampling than the DFT results. All adsorptions computed here indicate a negative Mulliken partial charge on the hydrogen ion, ranging from −0.10 to −0.30, indicative of surface-adsorbate charge transfer, similar to previous results. 4–6 Calculation of the changes in the surface work function from adsorption can be used to describe the formation of surface dipoles resulting from this charge transfer. 57 However, this is beyond the scope of our current effort. In order to examine the dynamic properties of our model, we have computed molecular dynamics trajectories at 298 K to determine the vibrational density of states for an H2 molecule adsorbed on either (100) and (111) surface (Figure 4). Ionic temperatures were controlled through use of Nosé-Hoover thermostat chains. 58,59 We have used a time step of 0.2 fs and Fermi-Dirac electron smearing 60 in order to avoid spurious forces due to negative electron occupation numbers from the Methfessel-Paxton smearing. The electronic temperature was set to 2000 K in order to aid in SCC optimization. Our calculations yielded approximately 0.75 ps/day for the (100) system and 1.25 ps/day for the (111) system, running with OpenMP threading on twelve 2.1 GHz Intel Xeon cores. We have computed total trajectories of up to 10 ps for each system. For the sake of comparison, we have also computed a DFT-MD trajectory for the (111) system with the Mermin functional and same electron temperature. In contrast, this yielded only 100 fs/day running on 288 Intel Xeon CPUs with the same clock speed, resulting in a trajectory of only 950 fs in total. This corresponds to an increase in computational efficiency of 102 − 103 for our DFTB model, similar to previous efforts. 25,27,28

13



Our results indicate that on the (100) surface, the H2 bond remains intact for at least some of the MD trajectory, as indicated by the peak at ∼4300 cm−1 . In contrast, the hydrogen molecule experiences rapid dissociation on the (111) surface, as indicated by the absence of this mode. In addition, the spectrum from the (111) surface indicates a slight blue-shifting of the H-Pu mode at ∼380 cm−1 , indicative of a slightly stronger H-Pu interaction due to the dissociation of the H-H bond. Our result for the (111) surface from DFTB compares well to the result from DFT, where we also see an absence of a H-H stretching mode and and H-Pu peak at ∼280 cm−1 . We note that this discrepancy in the H-Pu peak corresponds to an energy difference of only 0.015 eV, and is likely due in part to the poorly converged DFT-MD result. We compute an average distance of an H atom to the surface from DFTB of 3.62 Å for the (100) surface, and a value of 3.12 Å for the (111) surface. (The DFT result exhibited a long period oscillation between 2.9 and 4.7 Å.). The closer proximity of the molecule to the surface for the (111) system likely contributes to its enhanced dissociation rate. Next, we have computed minimum energy pathways (MEP) for surface dissociation and sub-surface diffusion using the climbing image nudged elastic band (NEB) approach. 61 All calculations were run within the Atomic Simulation Environment python interface, 62,63 with the image dependent pair potential method 64 used to generate an initial guess for each pathway. Here, our SCC convergence was held to 2.72 × 10−5 eV (10−6 au) and the NEB forces were converged to 10−5 eV. Calculations for MEPs with more than one energetic minimum were converged by performing additional calculations with the intermediate state as the initial/final NEB state. All transition states were verified through identification of a single imaginary frequency from a normal mode calculation, where the SCC convergence was held to 2.72 × 10−6 (10−7 au), and with hydrogen ion displacements of 10−4 Å. Analysis of the dissociation MEP on the (100) surface (Figure 5(a)) indicates a rate limiting step with an activation energy of 0.49 eV where the H2 molecule experiences a rocking motion before settling into a first minimum that is at an energy of −0.94 eV relative to the initial state.

14


Page 14 of 33

Page 15 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Here, one H ion is adsorbed into the (100) surface center site, and the second hydrogen ion adsorbed in a bridge site (i.e., between surface Pu atoms, approximately 1.20 Å above). The hydrogen ions then experience an additional rocking motion with an activation energy of 0.15 eV, before settling into the final state where both H ions are embedded within the center sites with a total reaction energy of ∆Erxn = −1.20 eV. These values are similar to our DFT energy grid calculations, where we determined a reaction energy of −0.85 eV for the first step, and total reaction energy of −1.22 eV. For dissociation on the (111) surface (Figure 5(b)), optimization of the initial state yielded a geometry with the H-H bond axis oriented 24.5◦ to the surface normal, a center of mass distance from the surface of 3.08 Å, and an adsorption energy of −0.05 eV, compared to DFT values of 31.1◦ , 3.34 Å, and −0.16 eV, respectively. We compute a small dissociation barrier of only 0.06 eV, which is commensurate with our MD results. In this case, the energy of the maximum in the MEP is approximately equal to that of the initial state, and a normal mode analysis yielded more than one imaginary frequency. We compute a total reaction energy of −1.30 eV, with a final state where the hydrogen ions are ∼0.80 Å above the surface in the hcp and fcc-type hollow sites. These results are similar to our DFT calculations, where we computed a total reaction energy of −0.95 eV. The difference between the calculations is likely due to the improved convergence criteria used in these calculations, as well as the approximate nature of the energy grid calculation in the DFT result, where hydrogen ions were placed in arbitrary positions on the Pu surface, rather than computing the NEB pathway. Finally, we have computed the MEPs for hydrogen ion subsurface diffusion on both δPu surfaces, which correspond to essential first steps in plutonium hydride formation. On the (100) surface (Figure 6(a)), we observe an activation energy of 0.92 eV for the initial reaction, which leads to an intermediate state that is at an energy of −0.94 eV relative to the initial state. The relatively low energy for the intermediate configuration is a result of the tetrahedral coordination of the hydrogen ion by four Pu atoms, where the H-Pu bond

15



distances are all approximately 2.0 Å. This hydrogen ion then experiences a relatively high activation energy of 2.22 eV before diffusing into the first sub-surface layer. The MEP is asymmetric due to the surface relaxation of the first two plutonium layers in our calculation. The final state of this MEP yields the H ion embedded within the sub-surface Pu layer, and is roughly isoenergetic with the initial state of the MEP. In contrast, on the (111) surface (Figure 6(b)) the hydrogen ion then experiences a relatively low activation energy of 0.66 eV, with a transition state where the H ion is coincident with the surface. The hydrogen then diffuses to an energetic minimum roughly equidistant from the top two Pu layers that is approximately isoenergetic with the initial state. It is likely that the MEP for diffusion through the next Pu subsurface layer will have a similar morphology to this part, with a slightly larger activation energy due to the aforementioned surface relaxation. These types of results can be used to determine ab initio kinetic Monte Carlo models for surface corrosion (e.g., Refs. 65–67 ), where libraries of possible reactions can easily be built due to the periodic nature of the crystalline slab and the relatively simplicity of hydrogen reactivity for this system.

4

Conclusions

In this work, we have created a DFTB model for hydrogen reacting with δ-Pu surfaces by parameterizing both the Hamiltonian and overlap matrices for the Pu-Pu and H-Pu interactions, and creating a many repulsive energy through use of the ChIMES Chebyshev polynomial formalism. Our new interaction potentials yield accurate results for the bulk δ-Pu electronic density of states, as well as the approximate H2 dissociation potential energy surface computed previously. We have computed vibrational spectra for H2 adsorbed on both the (100) and (111) surfaces, as well as minimum energy pathways for the dissociation reaction on each surface, where we observe somewhat more complex reaction paths than our previous results. We also have determined minimum energy pathways for sub-surface

16


Page 16 of 33

Page 17 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


diffusion of the hydrogen ion, where we find that the (100) system is possibly less prone to plutonium hydriding due to the presence of a stable, tetrahedrally coordinated inter-layer intermediate. The ChIMES interaction model allows us to include multi-center empirical DFTB repulsive terms in a straightforward way and with arbitrary complexity, to yield models that retain the majority of accuracy of DFT calculations with several orders of magnitude improvement in computational efficiency. These models can then be used to generate essential inputs for macroscopic models, such as kinetic data, which can then be used to simulate material aging and degradation (amongst other properties) with a high degree of accuracy.

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 borax and quartz massively parallel computers. The project 16-LW-020 was funded by the Laboratory Directed Research and Development Program at LLNL with N. G. as principal investigator. ChIMES potential parameters are available upon request.

17



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. Journal of Nuclear Materials 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 hydrogen-vacancy 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 2005, 72, 085101. (5) Huda, M. N.; Ray, A. K. A density functional study of atomic hydrogen adsorption on plutonium layers. Physica B 2004, 352, 5–17. (6) Huda, M. N.; Ray, A. K. An ab initio study of H2 interaction with the Pu (100) surface. Physica B 2005, 366, 95–109. (7) Y.Yang,; Zhang, P. Hydriding and dehydriding energies of PuHx from ab initio calculations. Phys. Lett. A 2015, 379, 1649–1653. (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. International Journal of Hydrogen Energy 2017, 42, 30727–30737. (9) Sun, B.; Liu, H.; H.Song,; 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. 18


Page 18 of 33

Page 19 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10) Söderlind, P. Ambient pressure phase diagram of plutonium: A unified theory for α-Pu and δ-Pu. Europhys. Lett. 2001, 55, 525–531. (11) Söderlind, P.; Landa, A.; Sadigh, B. Density-functional investigation of magnetism in δ-Pu. Phys. Rev. B 2002, 66, 205109. (12) Söderlind, P.; Sadigh, B. Density-Functional Calculations of α, β, γ, δ, δ 0 , and Plutonium. Phys. Rev. Lett. 2004, 92, 185702. (13) Söderlind, P. Quantifying the importance of orbital over spin correlations in δ-Pu within density-functional theory. Phys. Rev. B 2008, 77, 085101. (14) Liu, T.; Cai, T.; Gao, T.; Li, G. The electronic and structural properties of δ-Pu and PuO from the LSDA (GGA)+U method. Physica B 2010, 405, 3717–3721. (15) Söderlind, P.; Zhou, F.; Landa, A.; Klepeis, J. E. Phonon and magnetic structure in δ-plutonium from density functional theory. Scientific Reports 2015, 5, 15958. (16) 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. (17) Hernandez, S.; Freibert, F.; Wills, J. Density functional theory study of defects in unalloyed δ-Pu. Scripta Materialia 2017, 134, 57–60. (18) Goldman, N.; Morales, M. A. A First-Principles Study of Hydrogen Diffusivity and Dissociation on δ-Pu (100) and (111) Surfaces. J. Phys. Chem. C 2017, 121, 17950– 17957. (19) 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. 19



(20) Zheng, G.; Witek, H. A.; Bobadova-Parvanova, P.; Irle, S.; Musaev, D. G.; Prabhakar, R.; Morokuma, K.; Lundberg, M.; Elstner, M.; Köhler, C.; Frauenheim, T. Parameter Calibration of Transition-Metal Elements for the Spin-Polarized SelfConsistent-Charge Density-Functional Tight-Binding (DFTB) Method: Sc, Ti, Fe, Co, and Ni. J. Chem. Theory Comput. 2007, 3, 1349. (21) Koskinen, P.; Mäkinen, V. Density-functional tight-binding for beginners. Comp. Mater. Sci. 2009, 47, 237–253. (22) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931. (23) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a sparse matrix-based implementation of the DFTB method. J. Phys. Chem. A 2007, 111, 5678–5684, http://www.dftbplus.info (5 February, 2017). (24) Goldman, N. Multi-center semi-empirical quantum models for carbon under extreme thermodynamic conditions. Chem. Phys. Lett. 2015, 622, 128–136. (25) Goldman, N.; Koziol, L.; Fried, L. E. Using force-matched potentials to improve the accuracy of density functional tight binding for reactive conditions. J. Chem. Theory Comput. 2015, 11, 4530–535. (26) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. K. Construction of tight-binding-like potentials on the basis of density functional theory: Application to carbon. Phys. Rev. B 1995, 51, 12947. (27) Goldman, N.; Fried, L. E. Extending the Density Functional Tight Binding Method to Carbon Under Extreme Conditions. J. Phys. Chem. C 2012, 116, 2198–2204.

20


Page 20 of 33

Page 21 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(28) Goldman, N.; Srinivasan, S. G.; Hamel, S.; Fried, L. E.; Gaus, M.; Elstner, M. Determination of a density functional tight binding model with an extended basis set and three-body repulsion for carbon under extreme pressures and temperatures. J. Phys. Chem. C 2013, 117, 7885 – 7894. (29) Srinivasan, S. G.; Goldman, N.; Tamblyn, I.; Hamel, S.; Gaus, M. Determination of a density functional tight binding model with an extended basis set and three-body repulsion for hydrogen under extreme thermodynamic conditions. J. Phys. Chem. A 2014, 118, 5520–5528. (30) Koziol, L.; Fried, L. E.; Goldman, N. Using Force Matching To Determine Reactive Force Fields for Water under Extreme Thermodynamic Conditions. J. Chem. Theory Comput. 2017, 13, 135–146. (31) 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. (32) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953–17979. (33) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmentedwave method. Phys. Rev. B 1999, 59, 1758–1775. (34) Perdew, J. P.; Burke, K.; Enzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (35) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B 1993, 47, 558–561. (36) Kresse, G.; Hafner, J. Ab initio molecular dynamics simulation of the liquid-metalamorphous-semiconductor transition in germanium. Phys. Rev. B 1994, 49, 14251– 14271. 21



(37) 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. (38) Methfessel, M.; Paxton, A. T. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B 1989, 40, 3616–3621. (39) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188–5192. (40) 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. (41) Chen, S. P. Correlation-induced anomalies and extreme sensitivity in fcc Pu. Philosophical Magazine 2009, 89, 1813–1822. (42) Wang, Y.; Sun, Y. First-principles thermodynamic calculations for δ-Pu and -Pu. J. Phys.: Condens. Matter 2000, 12, L311–316. (43) Lashley, J. C.; Lawson, A.; McQueeney, R. J.; Lander, G. H. Absence of magnetic moments in plutonium. Phys. Rev. B 2005, 72, 054416. (44) 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.; Ramos, M.; Trouw, F.; Zhu, J.-X.; Haule, K.; Kotliar, G.; Bauer, E. D. The valencefluctuating ground state of plutonium. Sci. Adv. 2015, 1, e1500188. (45) Kullgren, J.; Wolf, M. J.; Hermansson, K.; Köhler, C.; Aradi, B.; Fauenheim, T.; Broqvist, P. Self-Consistent-Charge Density-Functional Tight-Binding (SCC-DFTB) Parameters for Ceria in 0D to 3D. J. Phys. Chem. C 2017, 121, 4593–4607. (46) Averill, F. W.; Painter, G. S. Steepest-Descent Determination of Occupation Numbers

22


Page 22 of 33

Page 23 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and Energy Minimization in the Local-Density Approximation. Phys. Rev. B 1992, 46, 2498–2502. (47) Hourahine, B.; Sanna, S.; Aradi, B.; Köhler, C.; Niehaus, T.; Frauenheim, T. SelfInteraction and Strong Correlation in DFTB. J. Phys. Chem. A 2007, 111, 5671–5677. (48) Köhler, C.; Seifert, G.; Frauenheim, T. Density-functional based calculations for Fe(n), (n ≤ 32). Chem. Phys. 2005, 309, 23–31. (49) Köhler, C.; Frauenheim, T.; Hourahine, B.; Seifert, G.; Sternberg, M. Treatment of Collinear and Noncollinear Electron Spin within an Approximate Density Functional Based Method. J. Phys. Chem. A 2007, 111, 5622–5629. (50) Press, W.; Flannery, B. P.; Teukolsky, S. A.; Wetterling, W. T. Numerical Recipes; Cambridge University Press: Cambridge, 1989. (51) 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. (52) Braams, B. J.; Bowman, J. M. Permutationally Invariant Potential Energy Surfaces in High Dimensionality. Int. Rev. Phys. Chem. 2009, 28, 577–606. (53) 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. (54) Aryanpour, M.; van Duin, A. C. T.; Kubicki, J. D. Development of a Reactive Force Field for Iron-Oxyhydroxide Systems. J. Phys. Chem. A 2010, 114, 6298–6307. (55) Behler, J. Constructing High-Dimensional Neural Network Potentials: A Tutorial Review. International Journal of Quantum Chemistry 2015, 115, 1032–1050.

23



(56) Behler, J. Perspective: Machine learning potentials for atomistic simulations. J. Chem. Phys. 2016, 145, 170901. (57) 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 2007, 75, 195112. (58) Nosé, S. Molecular Physics 1984, 52, 255. (59) Hoover, W. G. Physical Review A 1985, 31, 1695. (60) Mermin, N. D. Thermal properties of the Inhomogenous Electron Gas. Phys. Rev. 1965, 137, 1441–1443. (61) 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. (62) Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dulak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Jensen, P. B.; Kermode, J.; Kitchin, J. R.; Kolsbjerg, E. L.; Kubal, J.; Kaasbjerg, K.; Lysgaard, S.; Maronsson, J. B.; Maxson, T.; Olsen, T.; Pastewka, L.; Peterson, A.; Rostgaard, C.; Schiotz, J.; Schütt, O.; Strange, M.; Thygesen, K. S.; Vegge, T.; Vilhelmsen, L.; Walter, M.; Zeng, Z.; Jacobsen, K. W. The atomic simulation environment – a Python library for working with atoms. Journal of Physics: Condensed Matter 2017, 29, 273002. (63) Bahn, S. R.; Jacobsen, K. W. An object-oriented scripting interface to a legacy electronic structure code. Comput. Sci. Eng. 2002, 4, 56–66. (64) Smidstrup, S.; Pederson, A.; Stokbro, K.; Jónsson, H. Improved initial guess for minimum energy path calculations. J. Chem. Phys. 2014, 140, 214106. 24


Page 24 of 33

Page 25 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(65) Reuter, K.; Scheffler, M. First-principles kinetic Monte Carlo simulations for heterogeneous catalysis: Application to the CO oxidation at RuO2 (110). Phys. Rev. B 2006, 73, 045433. (66) Xu, L.; Campbell, C. T.; Jónsson, H.; Henkelman, G. Kinetic Monte Carlo simulations of Pd deposition and island growth on MgO(100). Surf. Sci. 2007, 601, 3133. (67) Mei, D.; Rousseau, R.; Kathmann, S. M.; Glezakou, V.-A.; Engelhard, M. H.; Jiang, W.; Wang, C.; Gerber, M. A.; White, J. F.; Stevens, D. J. Ethanol synthesis from syngas over Rh-based/SiO2 catalysts: A combined experimental and theoretical modeling study. Journal of Catalysis 2010, 271, 325 – 342.

25



Page 26 of 33

Table 1: Comparison of Selected Hydrogen Atomic Adsorption Energies and Approximate Distance from the Surface, with Comparison to Spin Polarization (SP) and Spin-Orbit Coupling (SOC) DFT Results. Surface (100)

(111)

Site Top

Method Ref 18 (SP) Ref 18 (SOC) Ref. 5 (SP) DFTB

Energy (eV) −1.87 −1.91 −2.12 −2.15

Distance (Å) 2.14 2.04 2.03 2.08

Bridge

Ref 18 (SP) Ref 18 (SOC) Ref. 5 (SP) DFTB

−2.66 −2.68 −3.14 −3.05

1.49 1.36 1.56 1.21

Center


−3.02 −2.91 −3.47 −3.25

0.00 0.00 1.13 0.00

Top


−1.75 −2.02 −2.52 −2.35

2.03 2.09 2.01 2.00

Bridge


−2.46 −2.55 −3.13 −2.50

1.57 1.52 1.57 1.67

hcp-type


−2.51 −2.61 −3.45 −2.71

1.53 1.52 1.42 1.31

fcc-type


−2.49 −2.60 – −2.86

1.49 1.48 – 1.50

26


Page 27 of 33

0.12

DFTB+U DFT+U (PBE)

0.1

0.08

Intensity

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.06

0.04

0.02

0 −3

−2

−1

0

1

2

Energy (eV) Figure 1: Total electronic density of states computed for δ-Pu for our DFTB model with spin-polarization and +U correction (red curve), with comparison to DFT+U (black curve). In each case, the density of states has been re-centered so the Fermi energy equals zero.

27



top top

center center bridge bridge

(a) (110) surface (a) (100) surface (a) (110) surface top hcp-type top hcp-type fcc-type fcc-type

bridge bridge (b) (111) surface (b) (111) surface (b) (111) surface

Figure 2: Renderings of the bare δ-Pu (100) and (111) surfaces, with labels for commonly identified adsorption sites. Here, plutonium atoms residing in the surface layer are colored black, atoms in the first sub-surface layer are red, and those in all subsequent layers are gray.

28


Page 28 of 33

Page 29 of 33

2.5

1 Pu

Pu Center site (dissociated)

0.5

1.5 1

0 Bridge site (dissociated)

Bridge site (dissociated)

0.5 Pu

−0.5

Center site (dissociated)

0 Initial state (H2 in center site)

Energy (eV)

[010] axis (Å)

2

−1 −1.5

1

2 3 [001] axis (Å)

4

(a) H2 dissociation potential energy surface

2.5

0.06 0.04

2

0.02 0

1.5

−0.02 −0.04

1

−0.06

Energy (eV)

[010] axis (Å)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


−0.08 0.5

−0.1 −0.12

0

−0.14 1

2 3 [001] axis (Å)

4

(b) Energy difference map from DFT+U

Figure 3: Map of the energy grid of H2 dissociation energies on the δ-Pu (100) surface, with comparison to DFT+U. Our DFTB model yields an accurate representation of the hydrogen surface chemistry for this system.

29



0.14 0.12 0.1

Intensity

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 33

0.08 0.06 0.04 0.02 0 0

500

1000

1500

2000

2500 3000 −1

3500

4000

Frequency (cm )

Figure 4: Power spectrum computed from DFTB-MD simulation of H2 on the (100) [red curve] and (111) [black curve] surfaces. Results from a 1 ps DFT+U-MD simulation on the (111) surface are shown with the thin black line. Both DFT+U-MD and our DFTB+U model indicate that the hydrogen molecule exhibits rapid dissociation on the (111) surface, leading to an absence of H-H stretch modes in the computed spectrum.

30


ii Adsorption Adsorption above center above center site site

0.80 0.80

Energy Energy(eV) (eV)

iii iii Upper H H ion ion Upper in bridging in bridging site site

iv iv

vv Both H H ions ions Both embedded in embedded in surface surface

ii ii

0.40 0.40 0.00 0.00

ii ii

0.49 eV eV EEaa == 0.49

ii

-0.40 -0.40 -0.80 -0.80

iv iv

iii iii

0.15 eV eV EEaa == 0.15

-1.20 -1.20

vv

-1.60 -1.60 (a) (100) surface

ii Adsorption above above Adsorption hcp-hollow site hcp-hollow site

ii ii

0.40 0.40

ii ii

ii


0.00 0.00

iii iii Both H H ions ions above above Both the surface surface the

-0.40 -0.40 -0.80 -0.80 iii iii

-1.20 -1.20 -1.60 -1.60 (b) (111) surface

Figure 5: Minimum energy pathways for H2 dissociation. The hydrogen atom initially closest to the surface is colored red for the sake of clarity. For the (100) surface, the final state indicates both hydrogen atoms adsorbed within the plutonium surface. In contrast, the energetic barrier for dissociation on the (111) surface is significantly lower, and the dissociated hydrogen atoms remain slightly above the plutonium surface. 31

Energy Energy (eV) (eV)

i i Surface adsorption Surface adsorption

1.60 1.60 1.20 1.20 0.80 0.80 0.40 0.40 0.00 0.00 -0.40 -0.40 -0.80 -0.80 -1.20 -1.20

ii ii Stable Stable intermediate intermediate

iii iii Sub-surface Sub-surface implantation implantation

Between Pu layers Between Pu layers

Ea = 0.92 eV Ea = 0.92 eV

iii iii Ea = 2.22 eV Ea = 2.22 eV

i i

ii ii (a) (100) surface

i i Adsorption above Adsorption surfaceabove surface


0.80 0.80

iii iii Interlayer Interlayer intermediate intermediate

ii ii

0.40 0.40 0.00 0.00

ii ii

Ea = 0.66 eV Ea = 0.66 eV i i

iii iii

-0.40 -0.40 (b) (111) surface

Figure 6: Minimum energy pathways for H diffusion to sub-surface layers on each plutonium surface. In this case, the (100) surface experiences higher diffusion energy barriers due to the formation of a tetrahedrally coordinated intermediate, where the hydrogen ion is simultaneously within a bonding distance to four plutonium atoms. The asymmetric MEP for the (100) surface is due to relaxation of the first two Pu layers in initial geometry optimization. 32

H2 dissociation potential 2.5

0

Bridging

1

−0.5

0.5

Pu

0

−1

Energy (eV)

1.5

Bridging

iii

0.60 0.20

−1.5 1

ii

0.5

Energy (eV)

2

i

Pu

Energy (eV)

Center site

Pu [010] (Å) [010] axis axis (Å)

Minimum energy path 1

-0.20 -0.60 -1.00

2 3 4 [001] axis (Å)

-1.40

TOC Graphic

33

i

ii iii

Development of a Multi-center Density Functional Tight Binding Model

Recommend Documents