Early Stage of the Dehydrogenation of NaAlH4 by Ab Initio Rare Event

Aug 28, 2012 - ... UCD-EMSC, University College Dublin, Belfield, Dublin 4, Dublin, Ireland ... Exploring the free energy surface using ab initio mole...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Early Stage of the Dehydrogenation of NaAlH4 by Ab Initio Rare Event Simulations Fabio Sterpone,*,† Sara Bonella,‡ and Simone Meloni*,§,∥ †

Laboratoire de Biochimie Théorique, UPR9080 CNRS, Institut de Biologie Physico-Chimique, 13 rue Pierre et Marie Curie, 75005 Paris, France ‡ Dipartimento di Fisica and CNISM Unit 1, Università “La Sapienza”, P.le A, Moro 2, 00185, Rome, Italy § School of Physics, Room 302 UCD-EMSC, University College Dublin, Belfield, Dublin 4, Dublin, Ireland ABSTRACT: We investigated the chemical species formed in the early stage of the first dehydrogenation reaction of (undoped) NaAlH4 (NaAlH4 ⇌ 1/3Na3AlH6 + 2/3Al + H 2 (↑)). We found that the experimental barrier to dehydrogenation (120 kj/mol) is compatible with the Al−H bond breaking of an AlH4 unit. We observed the formation of AlH3, AlH5, NaH, and NaH2. We computed the free energy profiles for the process of formation of the two most frequent species, AlH5 and NaH. While the free energy barriers for creating the species are comparable, the sample containing NaH is thermodynamically much more stable than the one containing AlH5. We did not observe the formation of H2 nor of the other products of the complete reaction. We attribute the lack of formation of H2 at this early stage to the fact that the hydrogen released in the sample is negatively charged and cannot quickly oxidate in the absence of chemical species that can efficiently be reduced. This finding suggests a possible mechanism for the catalytic action of Ti in Ti-doped samples. Finally, we studied the concentration of the various AlHi species as a function of the distance from the surface and found that species with higher negative charge are formed far from the surface, while neutral species are formed preferably closer to the surface, in particular in the top layer.

1. INTRODUCTION

Na3AlH6 ⇌ 3NaH + Al +

The development of a safe and efficient storage technology for hydrogen is a key factor for its practical use as an energy medium in many sectors, including the automotive one.1 Hydrogen can be stored in three ways: (i) as a compressed gas, (ii) as a cryogenic liquid, or (iii) in solid state materials. The last option is more convenient as it is safer and the hydrogen concentration is typically higher than in the first two. Materials for efficient solid state hydrogen storage must satisfy targets for the hydrogen content (expressed as volume and weight percent with respect to the material), hydrogen uptake and release rate, and lifecycle (expressed as the number of charge−discharge cycles before the material becomes unusable). A review of the target properties and different types of materials for solid state hydrogen storage can be found in refs 2 and 3. One of the most studied materials for solid state hydrogen storage is the sodium aluminum tetrahydrate (NaAlH4). The interest for this material is due to the discovery that Ti-doped samples undergo reversible hydrogen release/uptake4 at operative temperatures for automotive applications. The complete dehydrogenation of NaAlH4 occurs through three steps NaAlH4 ⇌

1 2 Na3AlH6 + Al + H 2( ↑ ) 3 3 © 2012 American Chemical Society

NaH ⇌ Na +

1 H 2( ↑ ) 2

3 H 2( ↑ ) 2

(2)

(3)

but only the first two are relevant as the third occurs at a temperature which is too high for practical purposes. Due to its low volume and weight percent of hydrogen, NaAlH4 cannot be used for technological applications. Nevertheless, sodium alanate is still the subject of intense research since it is considered a prototypical material from which valuable information can be obtained. In particular, a significant effort has been spent on the experimental and theoretical investigation of the mechanism of dehydrogenation of undoped and doped samples. For example, restricting our attention to some of the theoretical works more closely related to this paper, Ozoliņiš and co-workers5 investigated the kinetics, in terms of activation energy, of the Al mass transport via AlH3 and NaH vacancies. Wood and Marzari6 studied the existence of an unknown phase that might be favored at operating temperatures (∼400 K) over the ordinary lowtemperature AlH4 phase. The formation of this phase might Received: February 28, 2012 Revised: August 28, 2012 Published: August 28, 2012

(1) 19636

dx.doi.org/10.1021/jp3019588 | J. Phys. Chem. C 2012, 116, 19636−19643

The Journal of Physical Chemistry C

Article

precede the first dehydrogenation reaction. Vegge7 addressed the energetics of H2 desorption and the equilibrium structure of sodium alanate’s nanoparticles doped with Ti. Miceli and coworkers studied the mechanism and energy barrier of H2 release from the AlH4(101) surface catalyzed by a Ti4 cluster.8 Finally, several authors have tried to identify the chemical species formed during the dehydrogenation reaction both experimentally9,10 and by simulation.7,11−13 In this paper, we focus on the early stage of the dehydrogenation mechanism of NaAlH4 in an undoped sample. By the early stage we refer to the elementary steps that might initiate the reaction in eq 1. Our goal, in particular, is to identify, without making specific a priori hypothesis, the chemical species that might be formed as precursors of the products of this reaction and obtain some indication of the relative free energy cost of their formation. We also aim at identifying the distribution of the chemical species formed during the dehydrogenation as a function of the distance from the surface. This information might help to develop catalysts tailored to act on the chemical species present in the region around them. Even though we are restricting our attention to the early stages of the reaction, the time scales of the process are too long for brute force simulations (the reaction and the steps we focus on are rare events). To address the questions above, we thus combine ab initio simulations with a powerful approach for studying rare events, temperature-accelerated molecular dynamics (TAMD).14 This method allows us to accelerate the dynamics of a system by introducing a suitable set of collective variables (CVs) and coupling them to a set of dynamical variables that evolve at high temperature. In principle, TAMD can be used to reconstruct the free energy with respect to the CVs. However, this is computationally too expensive when the set is too large, like in the present case (see Section 2.3). Nevertheless, in this case TAMD can still be used to identify possible reaction paths, and we will use it in this spirit. As discussed more in detailed in the following, our simulations show that the breaking of Al−H bonds to form AlH3 and H− may not be the slowest step of the process. Rather, the further dissociation of AlH3 and the formation of H2 could be slower. These results are consistent with the hypothesis that the effect of Ti is to favor the recombination of H, more specifically H−, to form H2. They also seem to indicate that some significant structural change is needed to allow the progress of the reaction, as suggested by Wood and Marzari.6 The paper is organized as follows. In Section 2 we describe the method and the sample used in our simulations; in Section 3 we present and discuss the results; and finally, in Section 4 we draw some conclusions.

Figure 1. Sketch of the computational sample. The four bottom layers are kept frozen during the simulations for reproducing the effect of the bulk present in the real system. The six topmost layers are subject to the TAMD biased dynamics. The intermediate (buffer) layers are evolved according to unbiased Nosé-Hoover chain dynamics.

presence of the bulk. This setup has been validated by computing (1) the distance between the top two layers at 0 K (i.e., upon geometry optimization) as a function of the total number of layers and of the number of fixed ones and (2) the dehydrogenation energy (ΔE = Eslab − (Eslab−H + EH)) of an AlH4− unit in the top layer as a function of the amount of vacuum and of the width of the surface in the x−y direction. We found that for a system composed of 12 layers, of which 4 are fixed, the top interlayer distance is well converged. We also found that the ΔE converges very slowly with the amount of vacuum. This convergency can be significantly accelerated by using a modified Poisson solver that screens the electrostatic interaction of the system with its periodic images along the z direction. In particular, we found that by using the Martyna and Tuckerman Poisson solver15 the ΔE is well converged with 10 Å of vacuum above the top layer. Finally, we found that for a layer composed of 4 AlH4−/Na+ ions the dehydrogenation energy is well converged. 2.2. Temperature-Accelerated Molecular Dynamics Method to Study Rare Events. As mentioned in the introduction, the rare events method that we use, TAMD,14 is based on the introduction of an extended set of dynamical variables consisting of the atomistic degrees of freedom and an

2. ATOMISTIC MODEL OF THE NaAlH4 SURFACE AND COMPUTATIONAL METHODS 2.1. Atomistic Model of the NaAlH4 Surface. Following Vegge,7 we study the dehydrogenation process at the {001} NaAlH4 surface. Our atomistic model of the surface is composed of a slab of 2 × 2 × 1.5 unit cells, corresponding to 24 AlH4− and 24 Na+ ions. The system is organized in 12 alternate layers containing four AlH4− or Na+ ions each (see Figure 1). The supercell is orthorhombic with Lx = Ly = 10.016 Å and Lz = 27 Å. With this setup, we have ∼10 Å of unoccupied space in the z direction to model the vacuum above the surface. Four layers at the bottom of the slab are kept fixed to mimic the 19637

dx.doi.org/10.1021/jp3019588 | J. Phys. Chem. C 2012, 116, 19636−19643

The Journal of Physical Chemistry C

Article

reaction. While the choice of appropriate collective variables is always a delicate point when studying rare events, we believe to have included the most relevant chemical species for the reaction under investigation. Over a relatively short simulation (few picoseconds), we have also tested an extended set of collective variables including the number of molecules AlH and AlH2. This extended set did not bring any significant difference in the obtained results. In the following, we provide the detailed definition of the collective variables used in the present work. First of all, note that the chemical species considered here can be identified by calculating the coordination number with respect to H of the atoms Al, Na, and H. For example, the chemical species AlH4 is characterized by an Al atom 4-fold coordinated with hydrogen atoms. Similar considerations apply to the molecules belonging to the families NaHj and Hk. As a consequence, the number of molecules of a given species is defined as

auxiliary set of dynamical variables z. The z are associated to a set of collective variables θ(R) (R are the atomic positions, and θ = {θα} is a vector with dimensions equal to the number of collective variables) introduced to monitor the chemical reaction at hand. The dynamical variables z are a possible realization of the collective variables θ(R). The dynamics of the extended system is governed by the following coupled equations mR̈ = −∇R Uk(R, z) + thermo(β) μz ̈ = −∇z Uk(R, z) + thermo(β ̅ )

(4)

where m is a diagonal matrix whose elements are the masses of the atoms and thermo(β) indicates coupling to a thermostat at β = (kBT)−1. The tensor μ is the inertia of the variables z and is tuned so as to guarantee adiabatic separation between R and z; thermo(β̅) indicates that the auxiliary variables are coupled to a thermostat at β̅ = (kBT̅ )−1. The potential Uk(R,z) = V(R) + k/ 2(θ(R) − z)2 is the sum of the physical potential plus a term coupling R and z. In ref 14, it has been proven that under the condition of adiabatic separation between R and z, in practice with R moving much faster than z, and with k → ∞, the z move on the free energy surface of the atomic system at the physical temperature T, F(z). However, the auxiliary variables’ thermal energy is high enough to overcome barriers of order ∼ β̅−1 = (kBT̅ ). Thanks to the coupling with the z, the atomic system can also explore the configuration space much more efficiently than by brute force MD. In fact, TAMD is usually employed, choosing T̅ ≫ T, to quickly explore the free energy landscape and then combined with a suitable interpolation method16 to reconstruct it. While this is a powerful technique, it becomes too expensive when, as in our case, many collective variables are necessary to describe the reaction. In this work, we will employ the method in a different way. Several TAMD trajectories, each with increasing value of T̅ , will be run, and the evolution of the system will be monitored both by direct inspection of the atomic configurations and by analyzing the behavior of the collective variables along these trajectories. As we raise the temperature of the auxiliary variables, the system can overcome increasingly high free energy barriers. When these barriers correspond to bond breaking and/or formation in the system, new chemical species appear. Analysis of the trajectories then provides us with information on the identity and the “order of appearence” of the species in the early stages of the reaction and with an upper bound, ΔF ∼ (kBT̅ ), for their free energy of formation. 2.3. Collective Variables to Study the Dehydrogenation Process. The objective of this study is to investigate the mechanism of the early stage of the dehydrogenation reaction in terms of the chemical species formed during this process. It is therefore quite natural to use a set of collective variables that counts the number of molecules of the different chemical species expected to appear during the dehydrogenation that are effectively present in the simulation box. In particular, we have selected the following three families of chemical species: AlHi, NaHj, and Hk, with i = 3, ..., 6, j = 0, ..., 2, and k = 1,2. With this choice, we can accelerate and monitor the formation/breaking of the AlHi comprised between the reactant (AlH4) and product (AlH6) of the complete dehydrogenation reaction, and the species (AlH3) that we expect should be formed first in the reaction. We also monitor and accelerate the formation/ breaking of hydrogen molecules, and NaH and NaH2 species, that might be relevant intermediates or transition states in the

M

θα(R) = ν

∑ δl ,C (R) m

m=1

(5)

where α denotes the chemical species; l is the coordination number of the atom characterizing the species α (see eq 6); Cm(R) is the coordination number of the atom m; and δi,j is the Kronecker delta symbol. The coefficient ν is equal to one for all species but H2, for which ν = 1/2 to avoid double counting. In principle, the sum in eq 5 runs over all the atoms characterizing the species α. However, to avoid spurious effects due to the fixed atoms at the bottom of the slab, we have limited the sum only to the atoms belonging to the six top layers. This creates a buffer zone in which atoms evolve according to unaccelerated constant temperature MD in between the other two regions (see Figure 1). The function Cm(R) in eq 5 is defined as L

Cm(R) =

∑ (1 − Θ(|R m − R l|; R cut)) l=1

(6)

where Θ(x;Rcut) is the Heaviside step function, i.e., the function that is equal to 0 if x < Rcut and equal to 1 if x ≥ Rcut. In eq 6, the sum runs over the hydrogens in the sample. Obviously, in the case of H2, the H atom for which the coordination number is computed is excluded from the sum. With the definition above, our collective variables are not immediately suitable for TAMD simulations: since they are nonanalytical functions of R, we cannot compute their gradient with respect to the atomic positions ▽Rθα(R) as required by eq 4. We have solved this problem by introducing a smooth approximation to θα(R). This approximation is achieved in two steps. First, we introduce a smooth and continuous definition of the coordination number of an atom by replacing the term 1 − Θ(x;Rcut) with 1 − (x/Rcut)λ/1 − (x/Rcut)γ, where λ and γ are tunable parameters controlling the smoothness of the function. Second, we replace the Kronecker delta δl,Cm(R) with the Gaussian function exp[−(Cm(R) − l)2/2σ2]/(2πσ)1/2, where σ is a parameter controlling the accuracy of the approximation. In the rest of the manuscript we will denote the collective variables with the symbols {N(AlHi)}i=3,6, {N(NaHj)}j=0,2, and {N(Hk)}k=1,2, where in between parentheses we indicate the molecular species to which the collective variable refers. 2.4. Ab Initio Molecular Dynamics Simulation Setup. TAMD simulations were performed by running Density Functional Theory (DFT) based MD simulations as 19638

dx.doi.org/10.1021/jp3019588 | J. Phys. Chem. C 2012, 116, 19636−19643

The Journal of Physical Chemistry C

Article

implemented in the Quickstep module of the CP2K package.17,18 For the exchange and correlation terms we used the BLYP functional.19,20 The Kohn−Sham orbitals were expanded on a DZVP Gaussian Basis set, and the valence electron−nuclei interactions were modeled by the Goedecker−Teter−Hutter pseudopotentials.21 The electron densities were expanded on an auxiliary planewave basis set with a 280 Ry cutoff. The Gaussian-to-planewave basis set mapping was performed using the multigrid method with four grids. We used two different thermostats for the physical and collective variables’ subsystems, Nosé−Hoover chain22 and Langevin, respectively. We have chosen this approach as it was easier to implement the Langevin thermostat for the z variables. This mixed thermostatting was already successfully adopted in a previous work.12 The equations of motion of the z variables were integrated according to the Vanden−Eijnden−Ciccotti second-order algorithm.23 As for the parameter governing the dynamics, the Nosé−Hoover chain thermostat “frequency” (see appendix B of ref 22) was set to 2 ps−1, the inertia parameter μ to values in the range 1500−2000 amu bohr2, depending on the chemical species, and the friction parameter in the Langevin dynamics to 4.2 × 10 ps−1. The atomic masses were kept at their physical value. Finally, the time step was set to 0.25 × 10−3 ps.

Figure 2. N(AlHi) vs time at three different values of T̅ : T̅ = 6 × 103 K (top), T̅ = 12 × 103 K (center), and T̅ = 25 × 103 K (bottom).

AlH4 units (see red curves in Figure 4A). Initially, the distance between the Al of the donor (Ald) and the acceptor (Ala) units goes from ∼4 Å, the equilibrium distance between them, down to ∼3 Å. At this distance the H pointing toward the acceptor unit starts to experience the field created by the acceptor unit, and in fact, the amplitude of the vibration around its equilibrium position increases (see the top panel of Figure 4A). The distance of this hydrogen with respect to the acceptor Al oscillates following the Ald−Ala distance. When the H-bridge between the two units is established (t ∼ 0.25 ps in the figure), the oscillations of all the distances are strongly reduced. When the Al−H−Al bridge breaks, the system jumps back to its original structure. By increasing T̅ to 12 × 103 K, new chemical species are formed. In the central panel of Figure 2 we notice that the concentration of AlH4 decreases together with the increase of the concentration of AlH3. This seems to indicate that the experimental (free) energy barrier of ∼120 kJ/mol24 is due to the Al−H bond dissociation. As for AlH5, this kind of unit is formed only in a couple of events, and its stability is low (see events at t ∼ 12.5 ps and t ∼ 25−27.5 ps). This means that only a small fraction of the H atoms released by AlH4 units are absorbed by other AlH4 units. The other H released upon the Al−H bond breaking are captured by Na atoms, as shown by the value of N(NaHj) collective variables reported in Figure 5 (central panel). Analyzing in detail the concentration of the NaHj species we find that the most abundant is NaH. However, also a small but not negligible concentration of NaH2 is present in the sample. The mechanism by which a H is donated to an AlH4 unit to form AlH5 involves the formation of an intermediate structure in which the hydrogen atom is shared between the donor and the acceptor Al. A typical event of this kind is shown in Figure 3. In terms of interatomic distances, the process proceeds very similarly to the pure sharing described for the T̅ = 6 × 103 K. The main difference is that the process is terminated by a rotation of the acceptor unit that breaks the bond between the H and the donor Al. However, at this T̅ , after a few (1−3) picoseconds the system gets back to its initial state in which the donor and acceptor Al atoms are tetracoordinated.

3. RESULTS AND DISCUSSION In undoped NaAlH4 samples, the activation energy for dehydrogenation is estimated to be ∼120 kJ/mol,24 corresponding to ∼15 000 K. Such a high activation energy can originate from (i) the breaking/formation of H−Al bonds in AlHi units or (ii) the nucleation of separated products phases (Al and Na3AlH6). The collective variables described before and the size of the sample do not allow us to investigate the second phenomenon. Our setup aims instead at exploring the decomposition/formation path of AlHi units. As mentioned previously, for the time scale of the dissociation process to be on the time scale of MD simulations, the kBT̅ of the collective variables (see eq 4) must be of the same order of the relevant free energy barriers. We therefore run TAMD simulations at several T̅ temperatures in the range between about one-half and twice the value of the barrier reported above. In particular, we run three simulations at T̅ = 6 × 103 K, 12 × 103 K, and 25 × 103 K. In all simulations the physical temperature is kept fixed at T = 400 K, i.e., close to the temperature at which dehydrogenation is observed in undoped NaAlH4 (353 K25). We start the analysis of our simulations by comparing the curves N(AlHi) and N(NaHi) vs time at various T̅ . These data are shown in Figures 2 and 5, respectively. For T̅ = 6 × 103 K, we observe three events, two between 0 and 5 ps and one at about 12.5 ps, in which the number of AlH4 and AlH5 changes. The other (minor) changes in the concentration of AlH4 are due to the smoothed definition of the collective variables described in Section 2.3. Let us focus on these three “reactive” events. They are characterized by a decrease of the concentration of AlH4 and an increase of concentration of AlH5. Since no formation of AlH3 is observed, we conclude that the underlying process is hydrogen sharing between two AlH4 units. At such a low T̅ , the Al−H bond between the hydrogen and the Al atoms of the donor unit can not be broken, and the system quickly goes back to its original structure. This is confirmed by the visual inspection of the atomic trajectory. Examining in more detail the H sharing process, we observed that it is triggered by the reduction of the distance between two 19639

dx.doi.org/10.1021/jp3019588 | J. Phys. Chem. C 2012, 116, 19636−19643

The Journal of Physical Chemistry C

Article

Figure 3. Snapshots along a hydrogen transfer event between a pair of AlH4. The hydrogen is initially shared among two AlH4 units (a−b), thus forming an intermediate Al2H8. Then, at higher T̅ (T̅ = 12−25 × 103 K), the hydrogen is eventually donated from one AlH4 molecule to another (c), forming the pair AlH3 and AlH5. At T̅ = 25 × 103 K this process proceeds further, bringing the release of one H atom in the sample (d). The graphical representation of the figure and of Figures 6 and 7 has been obtained using the HyperBall plugin.40

Figure 4. Timeline of the donor−acceptor Al distance d(Ald−Ala) (bottom), Al acceptor−H d(Ala−H) (middle), and Al donor−H d(Ald−H) (top) for a hydrogen sharing event (A) and a hydrogen transfer event (B). The distances are expressed in Å, and the time is taken at an arbitrary origin close to the event. In the left panel we report two sharing events sampled at T̅ = 6 × 103 K and (red) and T̅ = 25 × 103 K. In the right panel we report a hydrogen transfer event sampled at T̅ = 25 × 103 K.

Finally, at T̅ = 25 × 103 K the AlH3 molecules formed are rather stable (i.e., once formed they do not transform immediately back to AlH4). The hydrogen released in the sample is partly bound to other AlH4 units, giving rise to AlH5 and partly to Na (forming both NaH and NaH2), as observed at T̅ = 12 × 103 K. However, in this case the stability of the AlH5 and NaHj units is increased. At this T̅ , the donation of a H atom can proceed through two different mechanisms: (i) by ion sharing, as already described for the simulation at T̅ = 12 × 103 K or (ii) by bond breaking and formation. In the latter process, the two Al atoms are not bridged via the hydrogen sharing and, as illustrated in Figure 4B, they oscillate around their equilibrium distance all along the H transfer. The first step of this process consists of the breaking of the Al−H bond. Initially, the Ald−H distance has an oscillation of large amplitude (see t = 0−1 ps in Figure 4B) and then suddenly jumps at ∼2.5 Å. The H atom remains at this distance for about 5 ps during which it attempts a few times to recombine with the donor unit (at t ∼ 1.5 ps, t ∼ 2.3 ps, and t ∼ 4.25 ps). These events correspond to times at which the Ald−Ala distance gets slightly smaller than the equilibrium distance (see the bottom panel of Figure 4B). Eventually, the H atom is definitively transferred to the acceptor when the donor and the acceptor move far apart (cf. top and bottom panel of Figure 4B for t ≥ 6.25 ps).

Figure 5. N(NaHi) vs time at three different values of T̅ : T̅ = 6 × 103 K (top), T̅ = 12 × 103 K (center), and T̅ = 25 × 103 K (bottom).

Also the transfer of a H atom from an AlH4 unit to a Na atom occurs via the sharing of the hydrogen (see Figure 6). 19640

dx.doi.org/10.1021/jp3019588 | J. Phys. Chem. C 2012, 116, 19636−19643

The Journal of Physical Chemistry C

Article

Figure 6. Snapshots along a hydrogen transfer event between an AlH4 and a Na. The hydrogen is initially shared among an AlH4 unit and a Na atom (a−b). The hydrogen is eventually donated to the Na atom (c) thus forming the NaH species. In most of the processes observed in our simulation the NaH species is unstable, and the proton is either donated back to the original Al or (d) is passed forward to another AlHi unit.

Figure 7. H donation process from one AlH4 unit to a Na atom. The center of the Maximally Localized Wannier Function localized nearby the moving atom is also reported as a yellow cloud.

of the homolytic breaking of the Al−H bond. We addressed this problem by computing the MLWF of unrestricted Kohn− Sham orbitals (the calculation is performed in zero-spin polarization conditions). In this case, the Al−H bond is described by two MLWF centers, one for the spin-up and one for the spin-down electrons involved in the bond. In the heterolytic bond breaking cases mentioned above, the two MLWF centers will stay close to the H or Al, forming H− or H+, respectively. If one center stays close to the H and one close to the Al, then the bond breaking is homolytic, and atomic hydrogen (H• radical) is formed. The results of the above analysis performed a posteriori on selected reactive events generated by the TAMD trajectories at T̅ = 12 × 103 K and T̅ = 25 × 103 K reveal that the breaking of the Al−H bond is heterolytic, leading to the formation of H− (see Figure 7). The negative charge of the atomic hydrogen present in the sample could explain the lack of formation of molecular H2 in our simulations. In fact, two H− ions strongly repel each other and cannot meet to form a H2 molecule. Moreover, and even more importantly, the reaction 2 H− → H2 + 2e is an oxidation (half) reaction and to occur requires the complementary reduction (half) reaction. For this reaction to take place, it must be present in the sample and, in a convenient configuration for an efficient electron transfer, a unit that can get reduced (i.e., accept the electrons released in the oxidation). We checked whether such a state exists among those visited in our simulations. The only state in which two H− ions are close together is when they are bound to a Na+ ion, so forming a ̂ NaH2 unit. Indeed, due to the large HNaH angle (∼120°), even in this state the two hydrogens are not very close (∼3 Å). In few cases, we observed states in which the angle HNaH ∼̂ 90°, corresponding to a H−H distance of ∼2.3 Å, still much larger than the H−H bond length in H2 (∼0.75 Å). Moreover, sodium, given its redox potential, is not a good partner for the oxidation of H−. Indeed, in the reaction paths proposed in the literature, typically H2 is not formed by free H− present in the sample. Shortly we will discuss these reaction paths and how the chemical species we identified in our simulations fit into them. For the time being, we notice that the observation that free hydrogen atoms do not form H2 suggests a possible mechanism for the formation of H2 in Ti-doped samples. Ti3+,

As discussed above, the TAMD simulations at the two higher temperatures reveal that the most frequent events are the breaking of the Al−H bond of an AlH4 unit and the donation of the hydrogen atom to either another AlH4 unit or to a Na atom. In the following, we will denote the first process as DAlAl and the second as DAlNa. To further characterize these two (competing) processes we compare their free energy profile as a function of a suitable collective variable. Following Monteferrante et al.,13 we use the difference between the distance of the migrating hydrogen from the donor and the acceptor: θ(R) = RDH − RAH, where RDH and RAH denote the first and second distance, respectively. The free energy profiles have been obtained by thermodynamic integration of the mean force that, in turn, has been computed via restrained molecular dynamics, as explained in ref 13. In the integration, we used 15 points for the DAlAl process and 25 for DAlNa. In both cases, the restrained molecular dynamics for computing the mean force employed a 4 ps long trajectory. As a first remark, both final states are found to be metastable; i.e., they are local minima of the free energy. The difference of free energy between the initial and final state of the DAlAl process is ∼90 kJ/mol, while for the DAlNa process this difference is equal to ∼62 kJ/mol. This implies that at the operational temperature the probability of formation of NaH is overwhelmingly higher than that of AlH5. The free energy barrier, on the other hand, turned out to be ∼100 kJ/mol for both processes. This means that if the dehydrogenation process is rate controlled, which seems to be the case given the very high free energy barrier, these products are equiprobable. The fact that we did not observe the formation of H2 even though we used a collective variable, N(H2), to accelerate its formation, and in spite of the high concentration of atomic hydrogen in the sample at T̅ = 25 × 103 K, deserves a careful investigation. To address this issue, we studied first the charge of the H atoms present in the sample. This was done by analyzing the relative distance between the H and Al atoms with respect to the center of the corresponding Maximally Localized Wannier Function (MLWF)26,27 during a dehydrogenation event. If the center remains close to the H atom, this means that a H− is formed. On the contrary, if the center remains close to the Al, then a H+ atom is formed. Care must be taken in the analysis to allow for an adequate representation 19641

dx.doi.org/10.1021/jp3019588 | J. Phys. Chem. C 2012, 116, 19636−19643

The Journal of Physical Chemistry C

Article

with its 3+ charge, can more easily attract two H−, increasing the concentrations of states in which two hydrogens are (relatively) close to each other. Moreover, its redox potential (E0 = 1.25228) makes the reduction half-reaction Ti3+ + 2e → Ti+ compatible with the oxidation half-reaction of the H−. Since the role of Ti is just that of a catalyst, this reaction needs to be followed by one or more reactions reestablishing the original oxidation of the Ti. However, once H2 is formed the system cannot easily get back to its initial state, and the reaction can more easily proceed toward the formation of Na3AlH6 and Al. This explanation is consistent with the experimental observation that TiCl3-doped samples contain a higher concentration of H vacancy.29 We stress that the above mechanism is only a proposal based on the anionic nature of the H released in the sample. Its validation goes beyond the goal of this work and would require us to assess accurately the structure of Ti-doped AlH4 samples, as well as the oxidation state of Ti, two issues still debated.30−35 We now compare our results with previous investigations to discuss how they fit into existing proposals for the overall reaction mechanism. Dilts and Ashby36 propose that the formation of Na3AlH6 proceeds via a double addition reaction of NaH to NaAlH4. As suggested by Ozoliņs ̌ and co-workers,5 this reaction will produce a NaH vacancy into the NaAlH4 that can migrate to the NaAlH4/Al interface where it annihilates, forming Al and H2. Also, the formation of AlH3, resulting from the donation of a H− to an AlH4−, fits into the mechanism proposed by Ozoliņs ̌ and co-workers. In this case, an AlH3 unit formed at the NaAlH4/Al interface decomposes into Al and H2 leaving a defective NaAlH4 phase. The corresponding AlH3 vacancy-AlH5 − defect migrates to the NaAlH4 /Na 3 AlH 6 interface and annihilates incorporating the excess Na+ into the growing Na3AlH6 phase. These processes start from the production of NaH and AlH3, which are indeed the species we observe. However, as we have seen above, in the absence of a catalyst favoring the formation of H2, their lifetime is short, in particular shorter than the time scale of diffusion of the species suggested by Ozoliņs ̌ and co-workers. This might also explain why hydrogen release in undoped samples occurs at high temperature,37,38,29 close to or above the melting of the sample, where the NaH and AlH3 vacancy diffusivity is much higher than at 400 K, the typical temperature at which hydrogen release is observed in doped samples. The AlH3, AlH5, and NaH species observed in our simulations are also key elements in the reaction path proposed by Walters and Scogin.39 We now analyze the concentration of the various species as a function of their position along the slab. The objective of this analysis is to evaluate whether some of them, namely, AlH3, AlH4, and AlH5, are preferentially formed close to the surface rather than in the bulk. This analysis has been performed on the simulation at T̅ = 25 × 103 K. The probability density function to find an AlHi unit at a position z from the bottom of the simulation box (PAlHi(z)) is reported in Figure 8 together with the integral IAlHi(z) = ∫ 0zdsPAlHi(s) measuring the fraction of the AlHi unit up to the level z of the slab. The difference δIAlHx(n) = IAlHx(z″) − IAlHx(z′) between the values of IAlHx(z) before and after the n-th layer measures the amount of a given species in the layer. The value of δIAlH4(n) is almost constant in the three AlHi layers subject to TAMD, namely, δIAlH4(1) = 0.21, δIAlH4(2) = 0.21, and δIAlH4(3) = 0.22 (the layers are numbered starting from the bottom, the layer 1 being the first

Figure 8. Distribution of the AlHx species along the z axis perpendicular to the crystal surface.

subjected to TAMD). On the contrary, the value of δIAlH3(n) is significantly higher in the top layer (δIAlH3(3) = 0.14) than in the other two (δIAlH3(1) = 0.075 and δIAlH3(2) = 0.055). Finally, the amount of AlH5 is almost the same in the first and second layer (δIAlH5(1) = 0.026 and δIAlH5(2) = 0.024), while this species is not present in the top layer. The above results can be explained as follows. In all the layers but the top one, the AlHi units are placed within a “cage” formed by eight Na atoms. In the top layer this cage is incomplete in the direction of the surface. As a result, species carrying negative extra charges, such as AlH52−, are preferably formed in the bottom layer, where the higher number of Na+ screens better the charge unbalance. For the complementary reason, the neutral species (AlH3) is preferably formed in the top layer.

4. CONCLUSIONS In this work, we have investigated the chemical species formed in the early stage of the dehydrogenation reaction of (undoped) NaAlH4. We found that the experimental barrier of 120 kJ/mol is compatible with the Al−H bond breaking of an AlH4 unit. We identified two main species formed when these bonds are broken and characterized the corresponding free energy barriers. The formed species, NaH and AlH5, have been mentioned as likely precursors of the dehydrogenation in the literature. However, we did not observe either the formation of H2 or the other products of the first dehydrogenation reaction (NaAlH4 ⇌ 1/3Na3AlH6 + 2/3Al + H2(↑)). We attribute the lack of formation of H2 to the fact that the atomic hydrogen available in the sample is negatively charged and to the absence of chemical species that can quickly oxidate it. This observation brought us to propose a possible mechanism for the formation of molecular hydrogen in Ti-doped material. Finally, we investigated the concentration of the various AlHi species as a function of the distance from the surface and found that species with higher negative charge are formed far from the surface, while neutral species are formed preferably closer to the surface, in particular in the top layer.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]; [email protected]. Present Address ∥

CASPUR, via dei Tizii 6, 00185 Roma, Italy.

Notes

The authors declare no competing financial interest. 19642

dx.doi.org/10.1021/jp3019588 | J. Phys. Chem. C 2012, 116, 19636−19643

The Journal of Physical Chemistry C



Article

(29) Singh, S.; Eijt, S.; Huot, J.; Kockelmann, W.; Wagemaker, M.; Mulder, F. Acta Mater. 2007, 55, 5549−5557. (30) Chaudhuri, S.; Muckerman, J. T. J. Phys. Chem. B 2005, 109, 6952−6957. (31) Thomas, G.; Gross, K.; Yang, N.; Jensen, C. J. Alloys Compd. 2002, 330−332, 702. (32) Sun, D.; Kiyobayashi, T.; Takeshita, H. T.; Kuriyama, N.; Jensen, C. M. J. Alloys Compd. 2002, 337, L8−L11. (33) Iñ́ iguez, J.; Yildirim, T.; Udovic, T.; Sulic, M.; Jensen, C. Phys. Rev. B 2004, 70, 060101. (34) Brinks, H. W.; Jensen, C. M.; Srinivasan, S. S.; Hauback, B. C.; Blanchard, D.; Murphy, K. J. Alloys Compd. 2004, 376, 215−221. (35) Graetz, J.; Reilly, J. J.; Johnson, J.; Ignatov, a. Y.; Tyson, T. a. Appl. Phys. Lett. 2004, 85, 500−502. (36) Dilts, J. A.; Ashby, E. C. Inorg. Chem. 1972, 11, 1230−1236. (37) Jensen, C.; Gross, K. Appl. Phys. A: Mater. Sci. Process. 2001, 72, 213−219. (38) Kircher, O. J. Appl. Phys. 2004, 95, 7748−7753. (39) Walters, R. T.; Scogin, J. H. J. Alloys Compd. 2004, 379, 135− 142. (40) Chavent, M.; Vanel, A.; Tek, A.; Levy, B.; Robert, S.; Raffin, B.; Baaden, M. J. Comput. Chem. 2011, 32, 2924−2935.

ACKNOWLEDGMENTS The authors thank M. Chavent and M. Baaden for the provision of the not-yet publicly available HyperBall plugin used for the graphical representation of the atomistic sample (http://www.baaden.ibpc.fr/hballs/index.html) and thank the CASPUR Supercomputing center for the provision of computational facilities. F.S. acknowledges CASPUR Consortium for financial support. S.M. acknowledges financial support from the European Community under the Marie Curie Intra-European Fellowship for Career Development Grant No. 255406. S.B. acknowledges financial support from the Istituto Italiano di Tecnologia under the SEED project grant No. 259 SIMBEDD - Advanced Computational Methods for Biophysics, Drug Design and Energy Research.



REFERENCES

(1) Schlapbach, L.; Züttel, A. Nature 2001, 414, 353−358. (2) Yang, J.; Sudik, A.; Wolverton, C.; Siegel, D. J. Chem. Soc. Rev. 2010, 39, 656−675. (3) Jena, P. J. Phys. Chem. Lett. 2011, 2, 206−211. (4) Bogdanovic, B.; Schwickardi, M. J. Alloys Compd. 1997, 1, 253− 254. (5) Gunaydin, H.; Houk, K. N.; Ozoliņs,̌ V. Proc. Natl. Acad. Sci. 2008, 105, 3673−3677. (6) Wood, B. C.; Marzari, N. Phys. Rev. Lett. 2009, 103, 185901. (7) Vegge, T. Phys. Chem. Chem. Phys. 2006, 8, 4853−4861. (8) Miceli, G.; Guzzo, M.; Cucinotta, C.; Bernasconi, M. J. Phys. Chem. C 2012, 116, 4311−4315. (9) Palumbo, O.; Cantelli, R.; Paolone, A.; Jensen, C. M.; Srinivasan, S. S. J. Phys. Chem. B 2005, 109, 1168−1173. (10) Palumbo, O.; Paolone, A.; Cantelli, R.; Jensen, C. M.; Sulic, M. J. Phys. Chem. B 2006, 110, 9105−9111. (11) Voss, J.; Shi, Q.; Jacobsen, H. S.; Zamponi, M.; Lefmann, K.; Vegge, T. J. Phys. Chem. B 2007, 111, 3886−3892. (12) Monteferrante, M.; Bonella, S.; Meloni, S.; Vanden-Eijnden, E.; Ciccotti, G. Sci. Modell. Simul. 2008, 15, 187−207. (13) Monteferrante, M.; Bonella, S.; Ciccotti, G. Phys. Chem. Chem. Phys. 2011, 13, 10546−10555. (14) Maragliano, L.; Vanden-Eijnden, E. Chem. Phys. Lett. 2006, 426, 168−175. (15) Martyna, G. J.; Tuckerman, M. E. J. Chem. Phys. 1999, 110, 2810−2822. (16) Maragliano, L.; Vanden-Eijnden, E. J. Chem. Phys. 2008, 128, 184110. (17) Lippert, G.; Hutter, J.; Parrinello, M. Mol. Phys. 1997, 92, 477− 488. (18) Lippert, G.; Hutter, J.; Parrinello, M. Theor. Chem. Acc. 1999, 103, 124−140. (19) Becke, A. Phys. Rev. A 1988, 38, 3098−3100. (20) Lee, C.; Yang, W. T.; Parr, R. G. Phys. Rev. B 1988, 37, 785− 789. (21) Goedecker, S.; Teter, M.; Hutter, J. Phys. Rev. B 1996, 54, 1703−1710. (22) Martyna, G. J.; Klein, M. L.; Tuckerman, M. E. J. Chem. Phys. 1992, 97, 2635−2643. (23) Vanden-Eijnden, E.; Ciccotti, G. Chem. Phys. Lett. 2006, 429, 310−316. (24) Sandrock, G.; Gross, K.; Thomas, G. J. Alloys Compd. 2002, 339, 299−308. (25) Yebka, B.; Nari, G.-A. Mater. Res. Soc. Symp. Proc. 2003, 801, BB4.4. (26) Marzari, N.; Vanderbilt, D. Phys. Rev. B 1997, 56, 12847−12865. (27) Berghold, G.; Mundy, C.; Romero, A.; Hutter, J.; Parrinello, M. Phys. Rev. B 2000, 61, 10040−10048. (28) Vany, P. Handbook of Chemistry: Electrochemical Series. In Handbook of Chemistry and Physics, 90th ed.; Chemical Rubbers Company: Boca Raton, FL, 2009. 19643

dx.doi.org/10.1021/jp3019588 | J. Phys. Chem. C 2012, 116, 19636−19643