Effects of solvation shells and cluster size on the reaction of aluminum clusters with water

the oxidation of aluminum nanoparticle: Multimillion-atom reactive molecular dynamics Reaction of aluminum clusters, Al n ( n = 16, 17 and 18), with liquid water is investigated using quantum molecular dynamics simulations, which show rapid production of hydrogen molecules assisted by proton transfer along a chain of hydrogen bonds (H-bonds) between water molecules, i.e . Grotthuss mechanism. The simulation results provide answers to two unsolved questions: (1) What is the role of a solvation shell formed by non-reacting H-bonds surrounding the H-bond chain; and (2) whether the high size-selectivity observed in gas-phase Al n -water reaction persists in liquid phase? First, the solvation shell is found to play a crucial role in facilitating proton transfer and hence H 2 production. Namely, it greatly modiﬁes the energy barrier, generally to much lower values ( < 0.1 eV). Second, we ﬁnd that H 2 production by Al n in liquid water does not depend strongly on the cluster size, in contrast to the existence of magic numbers in gas-phase reaction. This paper elucidates atom-istic mechanisms underlying these observations. Copyright 2011 Author(s). This article is distributed under a Creative Commons Attribution 3.0 Unported License .


I. INTRODUCTION
Reaction of water with metal to produce hydrogen gas has been widely studied, with aluminumwater reaction being most intensely investigated. [1][2][3][4][5][6][7] Although the reaction, 2Al + 6H 2 O → 2Al(OH) 3 + 3H 2 , is exothermic, formation of an aluminum-oxide layer on the aluminum surface prevents continuous reaction. 8 In order to overcome this bottleneck, continual removal of the oxide layer has been attempted using various promoters such as hydroxides, 5 oxides, 9 and salts. 10 Unfortunately, none of these techniques has achieved a sufficiently fast rate of H 2 production for commercialization. 8 Nanotechnology has opened new avenues toward solving this problem. For example, the surfacearea to volume ratio is drastically increased, 11 quantum effects on the atomic scale often lead to novel catalytic activities 12 and enhanced reaction rates, 1,7,[13][14][15][16] and surfaces can be tailored at the nanometer scale to prevent oxidation. 17 Metal clusters have been found to react chemically in unusual ways. A remarkable example is an Al superatom, i.e., a cluster consisting of a magic number of Al atoms. 18,19 Hydrogen molecules were produced through gas-phase reaction of Al superatoms with water molecules, where Al 12 and Al 17 showed particularly higher reactivity compared with other cluster sizes. 20,21 The size selectivity of superatoms for gas-phase reaction with water molecules has been explained through geometric factors. 20,21 First-principles calculations show the existence of sites with Lewis base character and those with Lewis acid character on superatom surfaces. In the proposed H 2 -production mechanism (hereafter referred to as R-mechanism), water molecules are first split into protons and hydroxyl anions, which respectively bond to Lewis base and acid sites. Subsequently, a hydrogen molecule is produced from two surface protons that are bonded to two proximate Lewis base sites, with the energy barrier of about 1 eV. Our previous quantum molecular dynamics (MD) study investigated H 2 production of the magic-number Al n (n = 12, 17) in liquid water. 22,23 It is found that proton transfer along a chain of hydrogen bonds (H-bonds) between water molecules, i.e. Grotthuss mechanism, [24][25][26][27][28] plays an important role not only in splitting water into surface hydrogen and hydroxyl groups, but also in recombining hydrogen atoms into hydrogen molecules. However, the H-bond chain is surrounded by non-reacting H-bonds that form a solvation shell, and the role of the solvation shell for H 2 production has not been addressed in the previous study. Thus, the central questions are: (1) What is the role of a solvation shell formed by non-reacting H-bonds surrounding the H-bond chain; and (2) whether the high size-selectivity observed in gas-phase Al n -water reaction persists in liquid phase?
To answer these questions, we investigate the reaction of superatoms, Al n (n = 16, 17, and 18), with liquid water using quantum MD simulations. Here, size-selectivity is studied by comparing highly reactive Al 17 with less reactive Al 16 and Al 18 , as observed in gas-phase reaction. 20,21 The observed fast reaction rates and calculated low reaction energy barriers are explained by the same reaction mechanisms as in our previous study. 22,23 In addition, we find that the formation and breaking up of branching H-bonds, which connect the H-bond chain with the solvation shell, are crucial to facilitate proton transfer. This explains why the energy barrier varies strongly upon the rearrangement of surrounding water molecules, generally resulting in extremely low values (< 0.1 eV) in the presence of branching H-bonds. Thermal deformation of superatoms also modifies the arrangement of Lewis base and acid sites, driving some of the surface hydrogen atoms to diffuse around neighbor aluminum sites and lowering the energy barrier. The Al 18 system undergoes further oxidation, from Al n -OH 2 to Al-OH-Al to Al-O-Al, revealing the mechanism of the formation of aluminum-oxide layers coating aluminum particles. Due to the thermal deformation of Al n , as well as the delocalized nature of the reaction assisted by H-bond chains, Al n -water reaction in liquid water is less sensitive to the local geometrical arrangement of Lewis acid and bases sites on the superatom surface, and accordingly the size-selectivity is less profound.

II. METHOD OF CALCULATION
Our quantum-mechanical molecular dynamics (MD) simulation is based on density function theory (DFT) within the generalized gradient approximation (GGA). 29 The projector-augmentedwave (PAW) method 30, 31 is used to calculate the electronic states, where the plane-wave cutoff energies are 30 and 250 Ry, respectively, for the electronic pseudo-wave functions and the pseudocharge density. A preconditioned conjugate-gradient method 32, 33 is used to minimize the energy functional iteratively. The Brillouin zone sampling is performed at point. The projector functions are generated for O (2s, 2p), H (1s), and Al (3s, 3p, 3d) states. Nosè-Hoover thermostat 34, 35 is used for MD simulation in the canonical ensemble, where the equations of motion are numerically integrated by an explicit reversible integrator 36 with a time step of 13 a.u. (∼0.314 fs).
With periodic boundary conditions, in a cell of 12.58×12.58×18.87 Å 3 , four systems are studied in our MD simulations. They are superatoms Al n (n = 16, 17, and 18) immersed in 84 H 2 O molecules, as well as a reference system consisting only of 96 H 2 O molecules. A snapshot of the Al 16 simulation cell is shown in Fig. 1(a). The MD simulations are performed at a temperature of 773 K to accelerate reactions, so that as many events as possible are sampled during 20 ps of simulation for each system. Once all reaction mechanisms in the MD simulations are identified, the corresponding reaction energy barriers are calculated with the nudged elastic band (NEB) method, 37 which is capable of finding the minimum-energy path between given initial and final states. With the help of transition state theory, 38 the reaction rates at room temperature are then estimated from the energy barriers.

III. RESULTS AND DISCUSSION
To study the hydrogen-production process, the time evolution of the numbers of H-H and Al-O bonds are plotted in Fig. 1(b). We consider two atoms to be bonded when their distance is less than a cutoff distance R c during a prescribed bond lifetime of 12 fs, where R c for Al-O, Al-H, and H-H are 2.4, 2.0, and 1.0 Å, respectively. The number of H-H bonds shows that, within 20 ps, three H 2 molecules are produced each in the Al 16 and Al 17 systems and that six H 2 molecules are produced in the Al 18 system. In contrast, no H 2 molecule is produced in the pure water system. The production of a larger number of H 2 molecules in the Al 18 system is attributed to the increasing number of Al-(OH)-Al bridging bonds, which will be discussed later in this paper. We do not observe any stable bridging bonds in Al 16 and Al 17 systems. After the first H 2 is produced in the Al 17 system, we also observe that one aluminum atom of the Al 17 -superatom is detached due to the high temperature. The detached Al atom forms H H > Al < O H O H 2 , which has two Al-H bonds. The three Al-H bonds in the final Al 17 system include two from the detached Al atom. Thus the final Al 17 superatom has only one Al-H bond corresponding to a surface H atom, which is the same number as in the Al 16 and Al 18 systems.
Before we investigate the H 2 production mechanisms in details, it is worth summarizing the mechanisms found in our previous study. 22 Figure 2 is a schematic of these mechanisms. An H-bond chain provides a channel for proton transfer (known as the Grotthuss mechanism) not only in the initial water splitting but also in H 2 production. The proton-donating side of the H-bond chain is a water molecule adsorbed at a Lewis-acid site (labeled A in  where Al and Al respectively denote aluminum atoms at Lewis base and acid sites. In Eq. (1), the reacting H-bond chain length l is defined as the number of oxygen atoms directly involved in the proton transfer. Reaction A → D is the water splitting process: Reaction A → B is proton exchange between Al n -bonded water and hydroxyl, where Al is another Lewis base site. In Eq. (3), the reaction could occur in both directions. Our previous study has revealed that the H-bond chain can substantially reduce the energy barrier of the reaction in Eq. (2) from 0.42 eV to 0.20 eV, 22 providing an example of H-bond related catalytic properties of water molecules. 41 To study the effect of H-bond chain length, we calculate the energy barrier for the H 2 -production reaction in Eq. (1) involving Al n (n = 16, 17, 18) with various chain length l = 1∼4. We start with a superatom in the minimum-energy state. The initial state of the reaction is prepared by first adding a water molecule to one of the strongest Lewis acid sites (i.e. the most negative site according to Mulliken population analysis) 13,42 of Al n to form Al-OH 2 , and adding a H atom to the nearby Lewis base site of Al n to form Al -H (see Fig. 3(a) for the case on Al 18 ). We then insert l-1 water molecules chained by the H-bonds to connect the Al-OH 2 and Al -H. The final state is prepared by proton transfer along the H-bond chain as shown in Fig. 3(a). Both initial and final states are relaxed to local minimum-energy configurations by the conjugate gradient method. Given the initial and final states, the minimum energy reaction path between them is obtained by NEB calculation. The difference between the maximum energy along the reaction path and the initial-state energy is the activation barrier of the reaction. Figure 3(b) shows the calculated energy barriers for Al n (n = 16, 17, 18) with various chain length l = 1∼4. The increasing rate of energy barrier versus length is about 0.20 eV for all three systems, which may be interpreted as the energy required for transferring a proton from a hydrogen donor to an acceptor. For l >1, we have obtained different energy barriers (within discrepancy ∼0.1 eV) even for the same superatom with the same H-bond chain length, depending on the initial arrangement of the water molecules along the H-bond chain.
Besides the proton transfer along the H-bond chain, the existence of massive H-bonds between the reacting H-bond chain and surrounding water molecules in liquid water (which we call the branching H-bonds) could play an important role in these reactions. To study how the H 2 -production mechanism is influenced by the branching H-bonds, Fig. 4 shows different stages of the H 2 production mechanism in the Al 16 system during MD simulation, for which the H-bond chain length l = 2. Figure 5 plots the distance R ij and bond-overlap population O ij between several atomic pairs i and j involved in this reaction. Here, the bond-overlap is calculated by population analysis 13,42 to clarify the change in the bonding properties. By expanding the electronic wave functions in an atomic-orbital basis set, we obtained the overlap population O ij between the i th and j th atoms as well as the gross charge Q i for the i th atom. The Al atom labeled Al1 in Figs. 4 and 5 is the proton-accepting Lewis base site, while Al2 is the proton-donating Lewis acid site. The reaction begins with the configuration in Fig. 4(a), where the O3· · ·H3-O1 H-bond has formed. Here, · · · denotes a hydrogen bond, and the H-bonds represented by red dashed lines in Fig. 4 are defined by the criteria: R oo < 3.5 Å and θ HOO < 30 • . The existence of those H-bonds has been confirmed by calculating the bond-overlap population O ij between the proton and the oxygen, as will be shown in the analyses below. The proton H4 transfers from O2 to O1 during (a)∼(b). As shown in Fig. 5(a), the bond-overlap between O3 and H3 increases, indicating the strengthening of the O3· · ·H3-O1 H-bond. At the same time, R O3-H3 decreases and R O1-H3 increases as shown in Fig. 5(b), indicating that H3 moves away from O1 towards O3. As H3 moves away from O1 and the H3-O1 bond weakens, it becomes for O1 to receive another hydrogen labeled H4. During (b)∼(c), breakage of the O3· · ·H3-O1 H-bond is indicated by the decrease of the bond-overlap population O O3-H3 in Fig. 5(a). Together with the decrease of O1-H3 and O1-H4 bond lengths, the breakage of the O3· · ·H3-O1 bond promotes breakage of the O1-H2 bond (R O1-H3 increases to ∼ 2 Å as shown in Fig. 5(b)) and pushes H2 into H1's binding region (R H1-H2 decreases to below 1 Å as shown in Fig. 5(b)). The O3· · ·H3-O1 bond breaks at 15.81 fs by way of water reorientation, which causes O3, H3, and O1 to misalign as shown in Fig. 4(c). This bond breakage is accompanied by the formation of another H-bond, O4· · ·H6-O3 in Fig. 4(c). Between (c) and (d), O1-H2 and H1-H2 bonds change until O1-H2 is broken and H1-H2 is bonded. From (d) to (e) is further structural change of the surface H 2 (H1-H2), and this H 2 molecule is released from the superatom surface during (e)∼(f). A notable observation during (c) and (f) is that the H 2 molecule is formed by the approach of H2 toward the Al1 atom while the H2-H1 bond is formed, instead of H1 being released to bond to H2 away from the superatom surface. We also observe that the O2-H5 bond within a surface-bonded hydroxyl group has higher kinetic energy (manifested by large-amplitude vibration of R O2-H5 ) than the O1-H3 and O1-H4 bonds in the free water molecule (see Fig. 5(b)). As demonstrated in the above example, branching H-bonds formed with surrounding water molecules (e.g. O3-H3 in Fig. 4) are actively involved in the proton transfer along the H-bond chain.
To quantify the effect of the surrounding water molecules (i.e. solvation shell) on the energy barrier for H 2 production, we consider the same H 2 -production mechanism as in Fig. 3(a) (with l = 1) by adding a surrounding water molecule; see H4-O2-H5 in Fig. 6(a). The extra water H4-O2-H5 serves as a proton donor to form a branching H-bond O2-H4· · ·O1. In order to study the dependence of the reaction energy barrier on the H4-O1 distance, we perform NEB calculation to obtain the H 2 -production reaction path, where R O1-H4 between H4 and O1 is constrained to a fixed value. The range of the constraint R O1-H4 is chosen in reference to the equilibrium O-H distance (1.75 Å) of H-bonds, which is the second peak position of the radial distribution function g(r) in our Al 18 MD simulation (see in Fig. 6(b)). (The first peak of g(r) signifies the covalent O-H bond length within water molecules.) By varying R O1-H4 = 1.75, 1.9, 2.0, and 2.1 Å, we obtain the energy barrier as a function of R O1-H4 as shown in Fig. 6(c). The energy barrier is reduced to below 0.02 eV (the blue dashed line) when R O1-H4 is 1.75∼2.0 Å, and the barrier increases to 0.08 eV for R O1-H4 = 2.1 Å. In the limit, R O1-H4 → ∞, we recover the case considered in Fig. 3 (i.e., in the absence of surrounding water molecules), where the energy barrier is 0.14 eV as indicated by the red dashed line in Fig. 6(c). This result indicates that branching H-bonds can greatly reduce the energy barrier of hydrogen production hence can increase the reaction rate.
To study the energy barrier with a larger solvation shell, we next perform NEB calculations using configurations from MD simulations, where the positions of the surrounding atoms (i.e., those outside the superatom and the reacting H-bond chain) are fixed. (A similar procedure of holding the surrounding water atoms to calculate the minimum energy path was used in Ref. 43.) Figures 7(a) and 7(b) respectively show the reactant configuration (corresponding to Fig. 4(a)) and the product configuration ( Fig. 4(f)) of the H 2 -production reaction in the MD simulation of the Al 16 system. For each configuration, we calculate the energy barrier of the same H 2 -production reaction involving Al 16 with the H-bond chain of length l = 2 using NEB, where the surrounding-atom positions are fixed. To prepare the initial state for the reactant configuration, we minimize the energy by relaxing the positions of the reacting H and O atoms (i.e. those in the reacting H-bond chain, which are colored in Fig. 7(a)) with all the other atoms fixed. To prepare the final state for the same reactant configuration, we move the reacting H and O atoms along a reaction path to form a product H 2 . The positions of the reacting H and O atoms are then relaxed to minimize the energy and obtain the final state for the reactant configuration. 44 Given the initial and final states, the reaction path with the surrounding atoms in the reactant configuration is calculated by NEB (see the red curve in Fig. 7(c)). The calculated final-state energy for the reactant configuration is 0.48 eV higher than the initial-state energy and is not stable. The initial and final states for the product configuration are prepared in a similar manner, and the calculated energy profile along the reaction path is shown by the red curve in Fig. 7(d). The reaction energy barrier for the product configuration is calculated to be 0.23 eV. The time difference between the reactant and product configurations in the MD simulation is ∼ 50 fs, which does not allow the superatom to deform noticeably. Therefore, the difference of Al 16 geometry between the reactant and product configurations is ignorable, and the difference in the energy barrier should arise from the different surrounding H-bond networks. The large energybarrier difference (from 0.48 to 0.23 eV) indicates the significance of concerted rearrangement of surrounding atoms to promote H 2 -production reaction. Since hydrogen molecules are not physically adsorbed to superatoms, holding non-reacting atoms fixed makes the H 2 confined closer to the superatom. Accordingly the final-state energy is higher than the initial-state energy (by 0.05 eV for the red curve in Fig. 7(d)).
To confirm the barrier reduction due to surrounding water in other systems, we also calculate H 2 production in the Al 18 system with the H-bond chain length l = 1 and 2. The corresponding energies along the reaction paths are respectively given by the blue and black curves in Figs. 7(c) (for the reactant configuration) and 7(d) (for the product configuration). For the reactant configurations, the energy barriers are typically less than 0.70 eV, while the product configurations are characterized by much lower energy barriers ( ∼ 0.05 eV). This again demonstrates the significant role of surrounding-atom configurations in lowering the energy barrier. The H 2 products in these two configurations are not tightly confined to the superatoms, and accordingly the calculated final-state energy is high (> 0.8 eV).
During our MD simulations, the Al 17 superatom splits into an Al 16 superatom plus one free Al atom, and the Al 18 superatom undergoes great deformation. These observations suggest that the geometry of the superatom plays an important role in the reactions. To investigate the effect of superatom deformation on H 2 production, Fig. 8 compares the energy barrier of H 2 production (see Eq. (4)) for Al n in the minimum-energy configuration and that for a deformed Al n taken from MD simulation with the H-bond chain length l = 1.
(4) Figure 8(a) is the initial and final states for the system involving the minimum-energy Al 18 . They are prepared in the same way as in Fig. 3(a). Figure 8(b) is the initial and final states for the system with a deformed Al 18 , of which atomic positions are taken from MD simulation and then O and H atom positions are relaxed to the local energy minimum. Figure 8(c) shows that the energy barrier is reduced from 0.14 eV (corresponding to Fig. 8(a)) to 0.004 eV (corresponding to Fig. 8(b)) due to the deformation of Al 18 . This is because the deformation of the superatom changes the effective charge of each site, in particular modifying the strength of the Lewis base site. This is reflected in the length of the Al-H bonds. For the minimum-energy Al 18 in Fig. 8(a), the Al-H bond length marked by the red circle is R(Al sym. -H) = 1.611 Å, which is shorter than R(Al def. -H) = 1.718 Å of the deformed Al 18 in Fig. 8(b). The same deformation effect is also observed for Al n (n = 16,17) as shown in Fig. 8(c). Besides the reduction of the energy barrier, deformation enhances the diffusion of surface H atoms on the superatom surface. Since a surface H atom bonded to the strongest base sites has the lowest energy, deformation drives the surface H atom to diffuse around several Lewis base sites, as the strength distribution of the Lewis base sites changes due to geometric deformation. This explains why we see migration of surface H atoms, reflected in spikes on the Al-H distance curve in Fig. 1(b). It is also worth noting that the Lewis-base character is quite sensitive to radicals (-OH etc.) bonded to the superatom. 20 Besides revealing plausible mechanisms of H 2 production, our MD simulation of the Al 18 system reveals the oxidization mechanism of Al atoms in water. The Al 18 -superatom undergoes higher geometric deformation, and hence the Al-(OH)-Al bridge bond is more favorable in an Al 18 system. Such bridge bonds seen in the Al 18 system are formed almost linearly depending on the oxidation level, i.e. the number of H 2 molecules produced, while on the contrary we did not observe stable bridge bonds in the Al 16 and Al 17 systems. Figure 9(a) through 9(c) shows the formation of bridge bond Al1-(O1H1)-Al2 in about 45 fs. The H-bond O1-H2· · ·O2 exerts force to drive O1 toward Al2, and the transfer of H2 from O1 to O2 occurs during the formation of the bridge bond. Subsequently, the O1-H1 bond in bridge bonds breaks to form Al1-O1-Al2 during Fig. 9(d) through 9(f). Al-O-Al bonds are basic building blocks of aluminum oxide, and the observed reaction from Al n (H 2 O) m to Al-(OH)-Al to Al-O-Al provides such an oxidation mechanism. This describes the formation of an aluminum-oxide layer, which covers aluminum particles to prevent the thermodynamically favorable reaction, 2Al + 6H 2 O → 2Al(OH) 3 + 3H 2 . According to the above oxidation mechanism, if the oxidation of Al-superatoms should be limited, the intermediate state, Al-(OH)-Al, should be avoided. The facts that Al 16 and Al 17 systems do not form Al-(OH)-Al bridge bonds and temperature-assisted geometric deformation facilitates such reactions provides a guideline for choosing the proper superatom size and reaction temperature for efficient H 2 production.

IV. SUMMARY
Our quantum MD simulation reveals H-bond assisted mechanisms for hydrogen production by the reaction of aluminum superatoms with water. Proton transfer along a chain of H-bonded water molecules is involved not only in the production of water but also in the splitting of water. Simulation results reveal the important role of solvation shell in accelerating the reactions. The NEB method is used to calculate various reaction paths, showing that H-bond chains as well as surrounding H-bonds in the solvation shell can substantially reduce the reaction energy barrier, thereby enhancing the reaction rate drastically. The character of Lewis acid and base sites is strongly dependent on the geometry of superatoms. However, thermal deformation of superatoms, as well as the delocalized nature of the H-bond-chain assisted reaction, makes the size-selectivity of the reaction less profound in liquid water compared with gas-phase reaction. Furthermore, the formation mechanism of aluminum-oxide layers coating aluminum in water is proposed.