A cluster density functional theory study of the interaction of the hydrogen storage system NaAIH4 with transition metal catalysts

(1)

A cluster density functional theory study of the interaction of the hydrogen storage system NaAIH4 with transition metal catalysts

Marashdeh, A.A.

Citation

Marashdeh, A. A. (2008, March 5). A cluster density functional theory study of the interaction of the hydrogen storage system NaAIH4 with transition metal catalysts.

Retrieved from https://hdl.handle.net/1887/12626

Version: Corrected Publisher’s Version

License: Licence agreement concerning inclusion of doctoral thesis in the Institutional Repository of the University of Leiden

Downloaded from: https://hdl.handle.net/1887/12626

Note: To cite this publication please use the final published version (if applicable).

(2)

Chapter 2 Methods and approximations

This chapter focuses on the theoretical methods employed in this thesis. In Section (2.1) the Born-Oppenheimer approximation is discussed. Section (2.2) gives a brief introduction to density functional theory (DFT). Section (2.3) deals with the geometry optimization method. In Section (2.4) the choice of a cluster model over a periodic model for representing the NaAlH₄ system is discussed.

2.1 The Born-Oppenheimer approximation: separating the nuclear and electronic motions

The Hamiltonian operator that describes both the nuclear and electronic motion is given by

H^∧tot(q_α,q_i)= K^∧ⁿ(q_α)+ K^∧^e(q_i)+V^∧^ee(q_i)+V^∧^ne(q_α,q_i)+V^∧ⁿⁿ(q_α) . (2.1)

In Eq. (2.1) Kn(q_α)

∧

and Ke(q_i)

∧

are the kinetic energy operators associated with the nuclei and electrons, respectively; Vnn(q_α)

∧

and V^∧ee(q_i) are the potential energy operators for the repulsions between the nuclei and the repulsions between the electrons, respectively; and V^∧ne(q_α,q_i) is the potential energy operator describing the attractive energy between the nuclei and the electrons. Furthermore, q_α and q_i symbolize the nuclear and electronic coordinates, respectively.

The total energy of the system can be calculated from the Schrödinger equation

H^∧totΨ(q_α,q_i)= EtotΨ(q_α,q_i) , (2.2) where Ψ is the total wave function of the system.

The solution of Eq. (2.2) can be simplified with an approximation that in many cases is highly accurate. Based on the fact that the nuclei are much heavier than the electrons, the

(3)

electrons move much faster than the nuclei and to a good approximation the nuclei can be considered as fixed in describing the motion of the electrons. Equation (2.1) can be written as

) (q_α K H Htot el n

∧

∧ = + , (2.3)

where the electronic Hamiltonian is given by

) ( ) , ( ) ( ) ( ) ,

(q_α q K q V q V q_α q V q_α

Hel i e i ee i ne i nn

∧

∧ = + + + . (2.4)

The Schrödinger equation for the electronic motion is then )

; ( ) ( )

; ( ) ,

(q_α q q q_α U q_α q q_α

H^∧el _i Ψ_el _i = Ψ_el _i . (2.5)

The energy U(q_α) in Eq. (2.5) is the electronic energy including internuclear repulsion, and represents the potential energy surface (PES) that determines the motion of the nuclei (see below). The electronic wave function depends parametrically on the nuclear configuration:

)

;

,_j(q_i qα el

el =Ψ

Ψ , (2.6)

where j labels the electronic state of the system.

Each time the electronic Schrödinger equation (2.5) is solved for a fixed set of nuclear coordinates, a “point” on the PES U_j(q_α) is provided. Whether calculated on the fly or represented through some fitting or interpolation scheme, this PES then determines how the positions of the nuclear coordinates will evolve through the time-dependent nuclear Schrödinger equation

) , ( ) ( ) ) (

,

( K q U q q t

t t

i nuc q n _nuc

α α

∂

∂ _Ψ

¸¹

¨ ·

©

§ +

Ψ = ^∧

. (2.7)

In Eq. 2.7 and throughout this chapter, atomic units have been used. The assumption that the electronic and nuclear motion can be separated, and that the nuclear dynamics takes place on a single PES representing a single electronic state, is known as the Born- Oppenheimer approximation [1]. Often it turns out that the dynamics of the nuclei can be

(4)

treated accurately through solving the classical equations of motion instead of the quantum mechanical equation in (2.7). However, in this thesis no dynamics of the nuclei is considered, the focus is on mapping out the most important and representative local minima on the PES.

2.2 Solving the electronic Schrödinger equation by Density Functional Theory In 1964 Pierre Hohenberg and Walter Kohn proved that for a nondegenerate ground state system of N interacting electrons in an external potential, the ground state energy of the system is uniquely determined by a functional of the ground state electron density [2]. In solving the electronic structure problem within the Born-Oppenheimer approximation [Eqs. (2.4) and (2.5)] the nuclei provide the external field [ V^∧ne(q_α,q_i) ]. By using the Hohenberg-Kohn variational principle [2], Kohn and Sham [3] found that the many–

electron problem [Eq. (2.5)] can be reformulated in a set of N single electron equations:

−1

2∇²+ν^(q)+ ρ^(q^'⁾ q− q^' dq^'

³

⁺^ν^xc^(q)

ª

¬

« «

º

¼

» » φi(q)=εiφi(q), (2.8)

where q and q^' in Eq. (2.8) are electron positions vectors; ρ(q^') is the electron density at q ;^' ν(q) is the external potential, which is related to V^∧ne in Eq. (2.1) through

V_ne = v(q)

³

ρ^(q)dr^{; and}^ν^xcis the exchange-correlation potential that is found from the functional derivative of the exchange-correlation energy E_xc[ρ(q)] with respect to the electron density:

v_xc(q)=δ^Exc[ρ^(q)]

δ ρ^(q) ^.^(2.9)

The electron density ρ(q) is given by ρ^(q)= φi(q)²

i

¦

N ^, (2.10)

where φ_i(q)^{is the i}^th one-electron Kohn-Sham wave function.

(5)

Equation (2.8) could be used to obtain the exact electronic ground state energyU₀(q_α), if the correct expression for the exchange-correlation energy functional E_xc[ρ(q)]^were known:

U₀(q_α)= εi i

¦

N ⁻¹₂

³

^ρ^(q)_q_{− q'}^ρ^(q^'⁾^dqdq^'^{+ E}^xc^[^ρ^]^{− v}

³

^xc^(q)^ρ^(q)dq^+Vⁿⁿ^(q^α^{) ,} ^(2.11)

where V_nn(q_α) is the nuclear-nuclear repulsion term described in Section 2.1.

The exact exchange-correlation functional described above has not yet been found, unfortunately, and only various approximations to it exist. One of these approximations is the local density approximation (LDA). This approximation employs the uniform electron gas formula for the exchange-correlation energy:

E_xc^LDA[ρ^(q)]=

³

ρ^(q)εxc(ρ^{(q)) dq}^.^(2.12) Here, ))ε_xc(ρ(q corresponds to the exchange-correlation energy per particle of a uniform electron gas of density ρ(q)^.

The LDA gives good results for quantities like bond lengths and vibrational frequencies of molecules, despite the rather severe approximation of treating the electron density as locally uniform. However, in determining molecular binding and adsorption energies, the LDA approximation is known to fail [4, 5, 6].

For this purpose, a more successful approximation for functionals is based on the generalized gradient approximation (GGA), where the exchange-correlation energy is written as a functional of both the electron density and its gradient:

E_xc^GGA[ρ^(q),∇ρ^(q)]= f (

³

ρ^(q),∇ρ^{(q)) dq}^.^(2.13) There are a number of different GGAs that give a considerable improvement over the LDA for energetics, such as the combination of the Becke correction [7] for the exchange energy and Perdew correction [8] for the correlation energy, or the gradient-corrected functional of Perdew et al. (PW91) [9].

(6)

In this thesis the binding energies of all the NaAlH4 clusters have been calculated using DFT [10, 11] as implemented in the ADF code [12]. In each case, a geometry optimization is performed (see Section 2.3), starting from a suitably chosen initial geometry based on the bulk crystal structure. The exchange-correlation energy is approximated at the GGA level using the PW91 functional [9]. The basis set used is of a triple zeta plus one polarization function (TZP) type. A frozen core of 1s on Al as well as Na was chosen, together with 1s2s2p for Ti and Sc, 1s2s2p3s3p4s3d for Zr and Pd, and 1s2s2p3s3p4s3d4p5s4d for Pt. The general accuracy parameter of ADF [12] was set to 4.0 based on a series of convergence tests. In many of the calculations a non-zero electronic temperature was applied to overcome problems with the SCF convergence.

However, it was ensured that the electronic ground state was reached by gradually cooling the electrons. The standard ADF fit sets (for the TZP basis sets) used to represent the deformation density were replaced by the fit sets corresponding to the quadruple zeta plus four polarization functions type basis sets. This was necessary since the standard fit sets were found to give inaccurate results. Results from tests showed that it was important to consider both spin restricted and unrestricted calculations. All calculations in this thesis have therefore been done both at the spin restricted and unrestricted levels. The spin unrestricted calculations were performed allowing one, two, three and four electrons to be unpaired. In the modeling of heavier transition metal catalysts (Zr, Pd, and Pt), scalar relativistic corrections were incorporated using the ZORA method [13].

In this thesis some periodic bulk calculations have been performed, those presented in Chapters 3 through 6 with ADF-BAND code [14] with the same basis and fit sets as used in the cluster calculations. The calculations presented in Chapter 7 were performed with the plane wave DFT codeVASP [15].

2.3 Finding (local) minima on the potential energy surface

A geometry optimization method can be used to find minima on a PES U₀(q_α), with these minimum energy structures representing stable (global minima) or metastable (local but not global minima) structures in which the system is most likely to be found.

(7)

The method used is based on a local second order Taylor expansion [16, 17, 18] of the PES U₀(q_α^l) around a given current structure ^q_α^l,

U₀(q_α^l+1)−U₀(q_α^l)=g^t⋅ Δq_α^l +1

2

( )

Δq_α^l ^t^{⋅ H ⋅ Δq}^α^l

1+ Δq

( )

_α^l ^t^{⋅ S ⋅ Δq}^α^l ^.^(2.14)

Here U₀(q_α^l⁺¹) is the estimated PES value at the next position ^q_α^l⁺¹, g^t the transpose of the gradient vector at q_α^l, Δq_α^l = q_α^l⁺¹− q_α^l, H the (approximate) Hessian (the second derivative matrix) at q_α^l, and S a matrix that allows for a scaling of the Hessian eigenvalues in order to follow all Hessian eigenvectors down-hill in energy towards a local minimum. Based on this Taylor expansion the best step along each eigenvector is determined to be

Δq_αl,i=− g_i λi−γi

, (2.15)

where g_i is the component of the gradient along Hessian eigenvector i, λ_i the corresponding Hessian eigenvalue, and γi is chosen to satisfy

λi−γi =1

2§ λi + λ2i + 4gi2

¹ ¸ . (2.16)

This method represents an iterative scheme in searching for a local minimum: First, an initial structure is chosen, the gradient calculated and an estimate made for the Hessian.

Then the Hessian is diagonalized and a step is taken along the eigenvectors according to Eqs. (2.15) and (2.16). Next, the gradient is calculated for the new geometry. If the largest component of the gradient and the previous step vector together with the change in energy are smaller than chosen thresholds, the geometry is considered to be converged.

The ADF default values have been used: 0.01 Hartree/Å for the components of the gradient, 0.01 Å for the components of the step vector, and 0.001 Hartree for the energy change from one structure to the next, respectively. If the geometry is not considered converged, the new Hessian is calculated based on the initial estimate and an added finite difference term built from the new and old gradient together with the step vector, and a step is taken along its eigenvectors. This process is repeated until convergence is reached.

(8)

2.4 A cluster approach to modeling NaAlH4

The interaction of an adsorbate with a metal substrate can be studied by modeling the substrate by a cluster with a limited number of substrate atoms [19, 20]. However, the interaction energy converges poorly with the size and shape of the cluster [19, 20]. This poor convergence behavior is associated with the metallic nature of the system. The problem can to a large degree be solved by using embedding techniques [21]. In modeling semiconductor materials, such as NaAlH4, the interaction energy converges much better with respect to cluster size and shape. The main reason is that the electrons are more localized.

Most of the theoretical work on modeling Ti-doped NaAlH4 has been using periodic calculations (see Chapter 1). These periodic calculations assume the system to be infinite, but from experiments it is known that the real NaAlH4 particles are nanometer-sized [22].

This is one of the reasons why the choice was made to work with a 23 NaAlH4 formula units (Z=23) semispherical cluster as a model for a nano-sized NaAlH₄ particle (see Chapter 3). The chosen cluster is about 2 nm large, which is small compared to a real NaAlH4 particle (typical particle sizes after ball milling are 150-200 nm [22]).

Nevertheless, the model cluster is structurally, electronically (density of states and band gap), and energetically (bonding energy per formula unit) close to bulk NaAlH₄ (Chapter 3). The cluster is chosen to have a large exposed (001) surface, because this has been shown to be the most stable surface of the different crystal faces [23]. In my opinion the cluster approach offers a number of advantages above the periodic approach. It is known that during dehydrogenation and hydrogenation large structural changes take place.

Modeling this properly might require very large periodic unit cells, making the calculations extremely demanding. By limiting the size of the unit cell to make the calculations feasible, one might introduce artifacts by allowing structural rearrangements to interact with their periodic images, thereby not obtaining the correct energetics.

Another point is that the real materials exhibit a range of surface facets, while a slab only has one. The edges and corners inherent in a cluster approach might therefore resemble the actual situation closer than when employing slabs. And, as already mentioned, the real NaAlH4 particles are nano-sized, as is my model cluster. It has been shown that when

(9)

the alanate particles are (significantly) smaller than 100 nm, NaAlH4 without Ti releases H2 at even lower temperatures than in Ti-catalysed NaAlH4 [24]. All this supports the view that a cluster model is worthwhile considering. Of course, I am fully aware that a cluster approach is not unproblematic, since the calculated properties might depend to some degree on both size and shape. However, since our common goal is to understand the role transition metals play in the dehydrogenation and hydrogenation of NaAlH4, I feel that all reasonable approaches should be explored. The cluster approach is certainly one of them, and one that could help provide important pieces of the puzzle.

2.5 References

[1] M. Born and R. Oppenheimer, Annalen der Physik, 84 457 (1927).

[2] P. Hohenberg and W. Kohn, Phys. Rev. 136 B864 (1964).

[3] W. Kohn and L. J. Sham, Phys. Rev. A 140 1133 (1965).

[4] J. A. White, D. M. Bird, M. C. Payne and I. Stich, Phys. Rev. Lett. 73 1404 (1994).

[5] G. Wiesenekker, G. J. Kroes, E. J. Baerends and R. C. Mowrey, J. Chem. Phys. 102 3873 (1995).

[6] B. Hammer, M. Scheffler, K. W. Jacobsen and J. K. Nørskov, Phys. Rev. Lett. 73 1400 (1994).

[7] A. D. Becke, Phys. Rev. A 38 3098 (1988).

[8] J. P. Perdew, Phys. Rev. B 33 8822 (1986).

[9] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B 46, 6671 (1992).

[10] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).

[11] W. Kohn and L. Sham, J. Phys. Rev. 140, A1133 (1965).

[12] G. T. Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders and T. Ziegler, J. Comp. Chem. 22, 931 (2001).

[13] E. van Lenthe, E. J. Baerends and J. G. Snijders, J. Chem. Phys. 101, 9783 (1994).

[14] G. T. Velde and E. J. Baerends, Phys. Rev. B 44, 7888 (1991).

[15] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).

[16] C. J. Cerjan and W. H. Miller, J. Chem. Phys. 75, 2800 (1981).

[17] A. Banerjee, N. Adams, J. Simons and R. Shepard, J. Phys. Chem. 89, 52 (1985).

(10)

[18] R. A. Olsen, G. J. Kroes, G. Henkelman, A. Arnaldsson and H. JÓnsson, J. Chem.

phys. 121 9776 (2004).

[19] M. A. Nygren and P. E. M. Siegbahn, J. Phys. Chem. 96, 7579 (1992).

[20] G. te Velde and E. J. Baerends, Chemical Physics 177 399 (1993).

[21] J. L. Whittem and H. Yang, Surf. Sci. Repts. 24, 55 (1996).

[22] H. W. Brinks, B. C. Hauback, S. S. Srinivasan and C. M. Jensen, J. Phys. Chem. B 109, 15780 (2005).

[23] T. J. Frankcombe and O. M. Løvvik, J. Phys. Chem. B 110, 622 (2006).

[24] C. P. Baldé, B. P. C. Hereijgers, J. H. Bitter and K. P. de Jong, Angew. Chem. Int.

Ed. 45, 3501 (2006).

(11)