Water Dissociation at the GaN(101̅0) Surface: Structure, Dynamics

C , 2012, 116 (27), pp 14382–14389. DOI: 10.1021/jp302793s. Publication Date (Web): June 21, 2012 ... C 2010 , 114 , 13695] have proposed a sequence...
0 downloads 11 Views 2MB Size
Article pubs.acs.org/JPCC

Water Dissociation at the GaN(101̅0) Surface: Structure, Dynamics and Surface Acidity Jue Wang,† Luana S. Pedroza,† Adrien Poissier,† and M. V. Fernández-Serra*,† †

Department of Physics and Astronomy, Stony Brook University, Stony Brook, New York 11794-3800, United States , and New York Center for Computational Science, Stony Brook University, Stony Brook, New York 11794-3800, United States



ABSTRACT: Pure GaN is known to show a very high photocatalytic water oxidation activity in the UV range. Recently Shen et al. [J. Phys. Chem. C 2010, 114, 13695] have proposed a sequence of intermediate steps for the water oxidation process at the GaN(1010̅ ) GaN/water interface. In this contribution we show that when spontaneous water dissociation occurs, the acidity of the surface can be accurately computed using first principles molecular dynamics simulations. The electronic structure analysis of the adsorbed water and hydroxyl groups shows large differences between GaN and the well studied photocatalyst TiO2. On the basis of our results we argue that the search for efficient photocatalytic materials needs to take into account the water dissociation activity of candidate material surfaces.



INTRODUCTION GaN:ZnO alloy semiconductors have been shown to be promising materials to serve as the photoanode in photocatalytical fuel cells.1 The microscopic details explaining origin of their high quantum efficiency (≈51% in the 420−440 nm range) to oxidize water in the presence of a sacrificial electron scavenger are unknown. However two physical reasons are highlighted as major contributors to the high photocatalytic oxidation activity of this alloy. The first one is its optimal band gap. As shown by Domen1,2 and subsequently studied by others, both theoretically3,4 and experimentally,5 the GaN:ZnO solid solution absorbs light in the visible range, while the optical activity of the two pure semiconductor counterparts only starts in the UV range. The second one is associated with the suitability of the conduction band edge (CBE) and valence band edge (VBE) potentials of the alloy for the overall water splitting reaction.6 In GaN the VBE, a band formed by nitrogen 3p states2,3 is located below the water oxidation potential,2 while the CBE, formed by hybridized Ga 4s and 4p orbitals is located above the water reduction potential. The activity of the solid solution as an overall water splitting material, despite the reduction in band gap, indicates that this band alignment with respect to the water redox potentials is maintained. There is a third and important question regarding the efficiency of a photocatalyst-material that has not been analyzed in detail to date. This is the importance of spontaneous water dissociation, H2O → H++OH−, at its surfaces. One of the most studied photocalysts, TiO2 is only active in the UV range of the solar spectrum. Many experimental and theoretical studies have shown that photoactive TiO2 surfaces do not dissociate water.7−9 On the other hand, the nonpolar GaN(101̅0) surface has been shown to spontaneously dissociate water.3,10 Even © 2012 American Chemical Society

more, very recently, it has been experimentally shown that the photocatalytic activity in this particular surface is very high.11 The question of how the efficiency of a photocatalyst correlates with the H2O dissociative potential of the surface, and how dissociated species interact with other water molecules at the interface, is yet to be addressed. In this contribution we analyze the dynamics of water dissociation at the GaN(101̅0). Our study represents, to our knowledge, the first thorough characterization of the changes in structure and dynamics of water interacting with a dissociating surface.



METHODS Ab Initio Molecular Dynamics. We have performed Born−Oppenheimer first principles molecular dynamics (FPMD) for a nonpolar GaN (101̅0) surface slab hydrated by water molecules in a supercell geometry as shown in Figure 1a. The slab supercell consists of 4 × 3 GaN unit cells repeated along the slab in plane directions. The short lattice constant obtained for the wurtzite GaN is a = 3.24 Å and the c/a ratio is 1.63 Å. Along the perpendicular direction to the slab surface (zdirection), periodic repetitions of a 5 layers slab (4 × 3 × 5 GaN units, 240 atoms) are separated by 16.2 Å of vacuum. This allows us to accommodate a total of 94 water molecules in our unit cell. To obtain an initial structure for the water layer, we cut a box of 94 water molecules with the appropriate geometry from a previous AIMD simulation of pure liquid water with 128 molecules.12 Therefore, our initial surface is fully nonhydroxylated at the beginning of the simulation. Received: March 23, 2012 Revised: June 13, 2012 Published: June 21, 2012 14382

dx.doi.org/10.1021/jp302793s | J. Phys. Chem. C 2012, 116, 14382−14389

The Journal of Physical Chemistry C

Article

3d states causes the cohesive properties such as lattice constant to contract about 2%, and has a larger effect on surface Ga−O bonds with around 4% reduction of bond length and 10% increase of bond energy. In our calculations, the pseudopotential for the bulk Ga atoms only includes the 4s2, 4p1 electrons as valence. However, for the double layer of surface Ga atoms, pseudopotentials with 3d10 as valence electrons are used. With these choices we are able to reproduce the absorption energies and geometries of the fully dissociated water monolayer on the surface described by Shen et al.3 with less than 0.3% error in the GaN lattice parameter and less than 7% error in the water adsorption energies. pKa Calculations: Theory. The computation of the logarithmic acid dissociation constant, pKa, is an important issue in electrochemistry.22 It can be obtained from the free energy for a proton diffusion process and different theoretical methods have been developed in the literature.22−26 In the present work, the pKa was derived from a method similar to that presented in ref 27 from the standard-state deprotonation free energy ΔG(0). The first step is the determination of the effective free-energy profile along the proton transfer coordinate. The choice of this coordinate is somewhat arbitrary, and should not influence the value of the free energy barrier. We choose to analyze the proton dynamics in the standard proton transfer coordinates ROO and δ. The coordinate ROO is the distance between the two O atoms (Oa and Ob), defined as the oxygen atoms that are involved in donating and receiving proton H. The δ is the displacement coordinate δ = ROaH − RObH of a given proton along the line joining the two corresponding oxygen atoms. The two-dimensional distribution function of the displacement coordinates, P(ROO,δ), is obtained as a normalized histogram evaluated from the corresponding FPMD coordinates. Integrating this distribution, we obtain the effective free energy profile ΔF(δ) as:

Figure 1. (a) GaN (101̅0) slab hydrated by 94 water molecules. (b) H+ and OH− bound to N and Ga atoms, respectively. N free surface atoms are also shown here. Liquid bulk water has been removed from the picture for more clarity. (c) OHA and OHD molecules chemisorbed on surface GaN.

The system was initially equilibrated at a temperature of 340 K using a Nose thermostat for 3 ps. After that, the FPMD simulation was continued within a microcanonical (NVE) ensemble using a Verlet integrator for a total of 10 ps, with a time step of 0.5 fs. This choice of temperature ensures the simulation of a liquid state of water with the basis set and functional of choice as shown in by Wang et al.13 The temperature conservation along the full trajectory is almost perfect, with a minor positive drift of 0.04 K/ps. The mean temperature along the 10 ps production run was 343 ± 10 K. It should be noted that the length of the simulation is not enough to ensure a complete sampling of the liquid potential energy surface, specially for the bulk region. However, as presented below, the solvation structure of the surface was equally reproduced on two different initial equilibrations which differ on the choice of pseudopotential for the Ga atoms (and also on the initial configuration of the liquid structure). In both cases we obtained the same overall dissociation results. The electronic structure and forces were obtained using density functional theory (DFT), within the generalized gradient approximation (GGA) and the PBE exchange and correlation functional.14 The Kohn−Sham wave functions were expressed as linear combinations of localized atomic orbitals as implemented in the siesta code.15,16 Core electrons were replaced by norm-conserving Troullier−Martins pseudopotentials17 in their fully nonlocal representation.18 A double-ζ polarized (DZP) basis set was used for both the semiconductor and the water molecules. As the covalent bond in a DZP basis is ∼0.5% stronger than in a larger and longer triple-ζ polarized basis,19 we expect that a DZP basis will slightly overestimate the water dissociation energy barrier, as compared to basisconverged results. As will be shown, this is not a problem, given the strong surface acidity obtained even with the use of this basis. Ga Pseudopotentials: The Role of 3d Electrons. The Ga pseudopotentials have a very important role on the structural and electronic properties of GaN, as has already been reported in refs 20,21. The Ga 3d electrons can be either treated as valence or as semicore electrons. These states strongly hybridize with the N 2s states near the bottom of the valence band and thus have a noticeable effect on the cohesive properties and energetic band structure of GaN. The neglect of

∫ dR OO P(R OO , R(δ))}

ΔF(δ) = −kbT ln{

(1)

The pKa is obtained from the deprotonation free energy ΔG(0). To do so, we map ΔF(δ) into ΔF(R), R being the ⎯⎯⎯→ proton displacement coordinate, or distanceOaH . The mapping is done by linear interpolation, by computing the mean ROaH at δ intervals of 0.02 Å width. With this, ΔG(0) depends on the free energy potential profile ΔF(R) as: ΔG(0) = −kBT ln{C0

∫R

R cut

dR A(R ) exp[−β ΔF(R )]}

min

(2) 2

where C0 indicates a 1.0 M concentration, A(R) = 4πR is a ⎯⎯⎯→ space phase factor determined by the integral of vector OaH ’s solid angle distribution times R2, and Rmin and Rcut are the distances (ROO and δ) delimiting the reaction. Finally, the proton affinity or pKa constant is computed as pKa = −log10 exp( −β ΔG(0))

(3)

where β is 1/kBT.



RESULTS AND DISCUSSION Surface Solvation. Due the periodic boundary condition used in the FPMD simulation, the top and the bottom GaN surfaces are wet, giving a total of 24 GaN units capable of

14383

dx.doi.org/10.1021/jp302793s | J. Phys. Chem. C 2012, 116, 14382−14389

The Journal of Physical Chemistry C

Article

region, and eventually return to a different surface Ga−OH, reforming a water molecule. In Figure 3, the average number of hydrogen bonds along the direction perpendicular to the GaN surface for: (a) surface

dissociating water. The surface water dissociation process is completed within less than 1.5 ps. During this time, surface Ga and N atom pairs catalyze the dissociation reaction until an equilibrium state is reached. In this way, we can distinguish a main type of dissociation events in this initial phase. The dissociation results in a OH− group bound to a Ga+ surface site and a H+ bound to the dangling bond of N− surface site. Three different types of surface bonds can be identified as the result of the reaction: surface Ga−OH, N−H bonds, and surface Ga− H2O bonds, formed by the molecular adsorption of water molecules to surface Ga sites. This can be observed in Figure 1b. As only ∼83% of the surface Ga−N pairs induce a full dissociation event, ∼17% of the surface N atoms maintain a dangling bond (N free in Figure 1b), free to receive a bulk water hydrogen bond (H bond). The initial dissociation process and the establishment of an equilibrium state is shown in Figure 2. This figure shows the

Figure 3. Average number of hydrogen bonds along the direction perpendicular to the GaN surface for (a) the surface OH groups, and (b) all water molecules. Curves display the total number of H bond (solid), the donor of H bond (dotted) and the acceptor of H bond (dashed).

OH− groups and (b) all water molecules are shown. We choose the following geometric criterion to define a hydrogen bond: oxygen−oxygen distance shorter than 3.5 Å and hydrogen bond angle larger than 140o. The Ga−OH bonds are strong enough to restrict hydroxyl groups to a narrow layer at the surface during MD simulations. At the two surfaces, the Ga-adsorbed water molecules are restricted to form donor type H bonds to bulk water because the lone pairs of the oxygen atoms are coordinated to surface Ga sites. The bulk water molecules maintain the nearly tetrahedral coordination of liquid water. Our simulation shows that about 80% of hydroxyl ions OH− are OH−−OH− dimers connected by H bonds. As a result we can distinguish the two components of the pair according to their donor/acceptor H bond character, as shown in Figure 1c. OHD is the hydroxyl group that donates its proton to its accepting neighbor OHA. The proton in OHA groups can form H bonds with bulk water molecules. However, this H bond with bulk water is formed only 75% of the time and is weaker than the H bond connecting the two hydroxyl ions. This effect is observed on the average O−H bond distance, (0.986 Å for OHA and 0.993 Å for OHD). Electronic Structure. In Figure 4, we present the electronic projected density of states (PDOS) onto selected atomic species. The plot is an average of 50 representative configurations along the full FPMD trajectory. The levels of the solvated GaN slab and those of the bulk GaN have been aligned with respect to the same vacuum level (at E = 0). However this alignment is subjected to fluctuations of the fermi level due both to solvent fluctuations and thermal vibrations of the slab. The bulk GaN Fermi level with respect to the vacuum level oscillated between −3.9 and −4.7 eV, these values are illustrated by the gray bar in the figure. The Fermi level is always located in the middle of the system’s band gap. It should be noted that the band gap of GaN as calculated with PBE is Eg ≈ 1.4 eV. Li et al.4 have shown that the use of GGA+U increases this gap by 1 eV. They also showed how this approximation mostly modifies the position of the conduction

Figure 2. Evolution of surface chemical bonds for the equilibration process (the first 1.5 ps) of FPMD simulation. The evolution of the total number of surface: (a) Ga−OH bonds, (b) Ga−H2O bonds, and (c) N−H bonds. The black curve represents the simulation performed considering the 3d states as valence electrons for the double layer of surface Ga atoms. The red one does not include them.

number of (a) Ga−OH, (b) N−H, and (c) surface Ga−H2O bonds as a function of time for two simulations (considering different Ga pseudopotentials). The geometrical distinction between Ga−OH and Ga−H2O groups is rather straightforward. We consider the H2O to be dissociated when a given OH distance is larger than 1.3 Å. The black curve represents the simulation done considering the 3d states as valence electrons for the double layer of surface Ga atoms. The red one does not include them (the same pseudopotential is used both for the bulk and for the surface Ga atoms). Although the dissociation rate can be affected by the choice of Ga pseudopotential, the full solvation process is well established within 1.5 ps, as shown in Figure 2. After this, the equilibrium phase still shows dissociation activity, observed in the fluctuations of Ga−OH and Ga−H2O bonds. These fluctuations correspond to two types of water dissociation, initiated from a surface Ga−H2O group, different to the process described above: (i) In some cases a surface Ga−H2O group will lose a proton to a neighboring Ga−OH group. This dissociation will leave the number of surface groups unaffected. However, it is possible to increase(decrease) the number of Ga−OH groups while decreasing(increasing) the number of surface-adsorbed molecules. (ii) This process occurs when a Ga−H2O group loses one proton to its neighbor hydrogen bonded bulk water molecule. This proton can now diffuse into the bulk water 14384

dx.doi.org/10.1021/jp302793s | J. Phys. Chem. C 2012, 116, 14382−14389

The Journal of Physical Chemistry C

Article

hopping diffusion in water. The details of the hopping and transport mechanism are still under debate both theoretically and experimentally.30−35 The main two proton states can be either H9O4+ (Eigen cation) or H5O2+ (Zundel cation). There are two plausible mechanisms for the proton transfer which are: (i) Eigen to Zundel to Eigen (EZE), as proposed on the basis of experimental NMR data33 and (ii) Zundel to Zundel (ZZ), as proposed on the basis of molecular dynamics simulations.30 Although proton diffusion is only possible in liquid water when excess protons exist, we have observed it during our FPMD simulation with the GaN(101̅0) surface serving as catalyst. From direct observation of FPMD simulations results, we have found that all proton diffusion processes are initiated from molecularly adsorbed water. As described above, H2O on the surface is bound to free Ga cation site by a “pseudo H bond” in which a Ga dangling bond attracts the lone pairs of O atom. This Ga−H2O group can be thought as a hydronium H3O+, with Ga replacing the excess proton. Water molecules in Ga−H2O groups can donate H bonds to neighbor OH− and to bulk water. Therefore, there are two possible pathways of proton diffusion along H bonds starting from Ga−H2O groups: (i) from surface H2O to neighbor surface OH−s (Diffs−s) and (ii) from surface H2O to bulk water (Diffs−b). After the hopping protons are transferred into bulk, we have also found proton diffusion through the standard Grotthuss mechanism: “Diffb−b”. Figure 5a−c shows these three types of proton diffusion

Figure 4. Projected density of states (PDOS), averaged over 50 different snapshots. Solid lines represent projection onto N orbitals: blue, surface N dangling bonds (N free); black, bulk N atoms (valence band edge). Non solid lines represent projections onto water related species: red dashed lines, surface hydroxyls (OH−); dashed green lines, H2O H bonded to dangling surface N; purple dot-dashed lines, bulk water; green dotted lines, Ga−H2O groups. The arrows at the top indicate the mean position of the top edge (as measured from zero to first maximum) of each corresponding projection. The gray bar indicates the fluctuations of the Fermi level of the system during the full FPMD. The Fermi energy is always located at the middle of the band gap.

band levels (Ga 3d states). Therefore, even if our band gap is largely underestimated, the following analysis of the PDOS should not be affected by this intrinsic PBE error on Eg, given that we limit our analysis to the valence band levels. According to the dissociation process already described, the water states can be separated into different groups: bulk water, water bound to surface Ga and water on top of the free N. We also present the PDOS for the hydroxyl ions at the surface. Although it is not shown in the figure (for clarity), the top edge of the OHA is always below the top edge of OHD, and their electric structure is noticeably different, indicating that they could be identified in XAS experimental studies of the surface. The levels of free N and bulk N are also shown in the Figure 4. The VBE of GaN corresponds to the top of the bulk N PDOS. In the water/TiO2 interface, the OH− levels can only lie above the VBE28 upon thermally activated solvent fluctuations which reduce the solvation shell of the ion. Otherwise the OH− ion is fully solvated, and its HOMO level sits below the VBE. In the present case, the OH− HOMO levels at the GaN(1010̅ ) are always ∼0.5 eV above the VBE of the bulk semiconductor. This indeed indicates that the OH− ions at the surface can receive the photogenerated holes. However, according to our results, the HOMO of the system will always sit at the non hydrogenated N atoms at the surface. These N dangling bonds usually are coordinated with a H2O molecule, which donates an H bond to the N. When this H bond is very strong, the HOMO can have a relatively large weight on this surface water molecule, as indicated by the large tail on the black curve (labeled “water top-free N” in Figure 4) which overlaps with the free N HOMO level. We can see that the water molecules bound to Ga appear to be the most stable in energy, as they can delocalize electrons through one of their lone pairs toward Ga, making a very strong accepting H bond, exactly as it occurs with the H bond acceptor in the water dimer.29 Proton Diffusion. In order to analyze the dynamics of proton diffusion, we should make connections to the Grotthuss mechanism, which is now a general notation for the proton-

Figure 5. Three types of proton diffusion processes are defined as (a) “Diffs−b”: proton diffusion between Ga−H2O and bulk water; (b) “Diffs−s”: proton diffusion between Ga−H2O and Ga−OH, and (c) “Diffb−b”: Grotthuss mechanism after hopping protons are transferred into bulk. The green atoms are oxygen atoms that are involved in proton transfer process, red atoms are oxygen atoms of neighbor water molecules, white are H, and brown are Ga atoms. Oa and Ob are defined in the text.

processes observed in the FPMD, where the green atoms are O atoms involved in proton transfer process, red atoms are O atoms of neighbor water molecules, H are white, and Ga atoms brown. We find that the frequencies of all types of proton transfer events after equilibration are rather high, indicating their free energy barriers should be comparable to room temperature thermal energy (kBT = 25 meV). Although the classical treatment of nuclei during proton transfer processes has been shown to overestimate effective free energy barrier by 50%,30,31 we can still compare the free energy barriers of different types of proton transfer processes in this study and compare those values with free energy barriers found for hopping excess protons in liquid water even though we are ignoring the zero point effects. As discussed in Figure 1, we need the effective free-energy profile (eq 1) for the proton diffusion process in order to obtain the pKa (eq 3). In Figure 6a−c, we show the twodimensional distribution function P(δ,ROO) (on the right side) and the effective free-energy profile (on the left side) for each 14385

dx.doi.org/10.1021/jp302793s | J. Phys. Chem. C 2012, 116, 14382−14389

The Journal of Physical Chemistry C

Article

Figure 6. Left panels are the effective free-energy profile along the proton transfer coordinate δ obtained from eq 1 for the proton diffusion processes: (a) “Diffs−s”, (b) “Diffs−b”, and (c) “Diffb−b” respectively. Right panels show the two-dimensional distribution function P(δ,ROO) normalized to unity and shown on the same scale.

duration of the simulation. The remaining four do not induce additional water dissociation and remain “free” of protons. In the Ga side, there are an average of four Ga forming Ga−H2O. These water molecules are the ones initiating all the deprotonation reactions described in the next sections. These protonation/deprotonation events are the source of the fluctuations in the total number of Ga−OH and Ga−H2O bonds seen in Figure 7. The number of events is large enough to allow us to compute the acidity of the surface. “Diffusion Surface to Surface” (Diffs−s). For the proton diffusion with initial Ga−H2O group from surface H2O to neighbor surface OH−, Oa is defined as the O atom of Ga−H2O and Ob is defined as the O atom of Ga−OH that receives the hopping proton H through hydrogen bond on the surface. Note that the definition of Oa and Ob can be switched, meaning Oa can be defined as O atom of Ga−OH and Ob can be defined as O atom of Ga−H2O. The reason is that once the proton is transferred, the original Ga−H2O becomes the new Ga−OH and vice versa. In this case, the distribution function P(δ,ROO) is characterized by two symmetric prominent peaks around (δ,ROO) ≈ (±0.5,2.5) Å as shown in Figure 6a. It also has significant non-negligible weight at δ = 0 Å, similar to what occurs in the H5O2+ (Zundel cation) in liquid water. For this diffusion process, the effective free-energy barrier of proton hopping between Ga−H2O and Ga−OH (symmetric process) is around 40 meV, which is lower than the free-energy barrier of “Diffs−b”.

process. The definition of Oa and Ob used in the displacement coordination δ is different for each of these processes and will be described below. In Figure 7 we show the number of N−H, Ga−OH and Ga− H2O bonds at the surface as a function of the simulation time. After 1 ps of simulation a total 20 dissociation events resulting in Ga−OH and N−H bonds occur. After this initial equilibration time, a dynamic equilibrium state settles. Among the 24 N atoms, 20 remain bound to a H atom for the full

Figure 7. Number of N−H, O−H, and H2O at the surface as a function of the simulation time. 14386

dx.doi.org/10.1021/jp302793s | J. Phys. Chem. C 2012, 116, 14382−14389

The Journal of Physical Chemistry C

Article

sections. The so-called “diffusion surface to surface” process (Figure 6a), does not contribute to the surface acidity, given that the proton hops through but never leaves the surface. From the previous results we see that the deprotonation free energy barrier for transferring a proton from a surface water to a bulk water (“Diffs−b”) at the GaN(101̅0)/aqueous interface is around 75 meV. Using eqs 1 and 3, we obtained pKa = 3.0 ± 0.1 from our FPMD simulations. This indicates that the GaN(101̅0) surface yields an highly acid aqueous medium. Even more, as observed from the deprotonation free energy barriers obtained in the three different processes described above, once the protons have left the surface, the barriers to diffuse within the bulk liquid region are smaller than those to return to the surface. We can apply the same procedure to obtain the pKa for this Diffb−b process. This process is nothing else than the protonation reaction of water H2O + H+ ⇌ H3O+, and the free energy profile is shown in Figure 6c. We obtain pKawater = 1.9 ± 0.1. The value is consistent with the concentration of protons we have within the bulk water region, which is approximately one proton in a bulk of 70 water molecules. For both pKa calculations, the largest error originates on the choice of Rcut in eq 2. R is the distance traveled by the proton from the reactant to the product. As our free energy profile is obtained from a dynamic equilibrium reaction, Rcut is chosen to be at the middle between the minimum in the product valley and the top of the reaction barrier. Shifting Rcut closer to the maximum increases the pKa up to 0.1 units. Analysis. There remains the question of how much these results would change if we would have started our FPMD simulations from a fully hydroxylated surface. In that case, if the results shown here represent indeed an equilibrium state, we would expect to see the deprotonation of some of the N−H surface groups and subsequent reformation of water molecules by protonation of some of the Ga−OH groups. The deprotonation events we observe in our simulation never result in the formation of N−H bonds; i.e., the unsaturated N atoms at the surface remain unprotonated. This occurs even if some of the deprotonation events initiate from water molecules adjacent to these free N sites. This indicates that the four free N surface atoms have a negligible proton affinity. We should expect this to be independent of the initial conditions. Indeed, the proton affinity we obtained is intrinsically related to both the charge of Ga and N atoms at the surface. To show it we need to relate our results to the physical meaning of surface acidity or pKa. One of the standard empirical methods to evaluate surface acidity is the modified version of the MUSIC model (multi site complexation model),36 which uses the Pauling bond valence to predict proton affinity constants or pKa values. This model was explicitly developed for oxide surfaces, but we can use it to estimate the charge of Ga atoms at the surface. According to the MUSIC model, pKa = −AQ, where A = 19.8 is an empirically fitted constant,36 and Q is the charge of the O atom at the surface undergoing the protonation. To obtain this charge, this model needs to estimate the coordination number of the protonated base,26 and justify the choice a posteriori by agreement with experimental data. If we apply this model in reverse to our calculated pKa, we obtain that the charge of the hydroxyl O atom being protonated is Q = −0.15 ± 0.005e−. This Q is the same that appears in the reaction H+ + (Ga−OH)Q ⇌ (Ga− OH2)1+Q. Therefore, Ga−O bonds at the surface contain an electron excess of at least 0.15 electrons.

“Diffusion Surface to Bulk” (Diffs−b). In this case, Oa is always defined as the O atom of H2O bound to surface Ga cation sites (Ga−H2O) and Ob is always defined as the O atom of the water that receives the hopping proton H through the hydrogen bond. The distribution for Diffs−b is characterized by two asymmetric peaks around (δ,ROO) ≈ (−0.55,2.6) Å and (0.7,2.6) Å in Figure 6b, showing the geometric character of OaHOb triplets. The P(δ,ROO) distribution of this type of proton diffusion also has non-negligible weight at δ ≈ 0 Å. This evidence of existence of centrosymmetric complexes (Ga− H2O···H2O) is similar to H5O2+ (Zundel cation) in liquid water. However, the distribution shows that protons are more likely to originate from Ga−H2O. Using eq 1, we can obtain that the effective free-energy barrier for proton transferring from Ga−H2O to bulk water is around 75 meV. The reverse path, from H3O+ to Ga−OH, has a free-energy barrier around 55 meV. “Diffusion Bulk to Bulk” (Diffb−b). Finally, we have also studied the Grotthuss mechanism after hopping protons are transferred into the bulk water region. It is worth recalling that Ga−H2O can be analogous to H3O+. Like H3O+, we should always find the hopping proton in the Ga−H2O. In our case, the transferring proton is defined as the one that has the smaller absolute value of δ, meaning that this proton is more likely to hop. Once the hopping proton is transferred into bulk, then it becomes an excess proton in liquid water. We use the same definition of Oa and Ob for “Diffb−b” as in ref 30, shown in Figure 6c. Oa is defined as the center O atom of the most probable H3O+, and there are three H atoms to form H bonds with neighbor molecules. Ob is defined as the O atom that has the smallest absolute value of δ, meaning this H atom is the most probable hopping proton. The free-energy barrier for “Diffb−b” is obtained to be around 20 meV, about the same as the value in a classical limit obtained in ref 30. As in excess protons in liquid water, the effective free energy barrier is also very small (comparable to kbT at room temperature). It has been shown that the rate of excess proton diffusion in liquid water is determined by thermally induced hydrogen-bond breaking in the second solvation shell. Surface Acidity: pKa. The acidity of a surface is a measure of its deprotonation potential, i.e. its ability to lose protons to bulk water. In our simulation, this is achieved through the socalled “diffusion surface to bulk” process (Figure 6b). The deprotonation process we observe here is somehow unusual. It initiates from a H2O molecule bound to a surface Ga atom, which is not hydroxylated. The reaction is H+ + Ga−OHQ ⇌ Ga−OH21+Q. The charge Q is unknown and depends on the actual charge the Ga atom at the surface. This unusual deprotonation process can be understood thinking on the surface as a catalyst for water autodisociation, releasing protons to the bulk liquid region and maintaining the corresponding hydroxyl ions adsorbed. To understand this better, we should remind that we did not start our FPMD from a fully hydroxylated surface, instead, we started from a pure GaN(101̅0), in contact with liquid water. The unit cell chosen is such that 12 GaN pairs are exposed to water at each side of the surface slab. As already explained, and shown in Figure 7, after 1 ps of simulation 20 dissociation events occur, leaving the surface with 4 GaN pairs free to interact with bulk water. The free Ga sites strongly bind water molecules. These molecules remain bound to the Ga for the full duration of the simulation, and they undergo two different types of deprotonation processes, as described in the previous 14387

dx.doi.org/10.1021/jp302793s | J. Phys. Chem. C 2012, 116, 14382−14389

The Journal of Physical Chemistry C

Article

Notes

In wurtzite semiconductors each atom contributes with 1/4 of its valency to each of its four bonds. In GaN, each Ga provides 3/4 e− and each N provides 1/4 e− per bond. The “electron counting rule”37 (ECR) for GaN states that surface dangling bonds are empty for Ga (Ga3/4) and full for N (N−3/4). We can make a simple calculation, taking into account our result and the ECR if we assume that the Ga−O bond has an excess of 0.15 e− and if we assume that this excess charge in the Ga−O bond (to which the OH− group contributes with 1 e−) is transferred from the surface N. In that case, the total charge available per N atoms is as follows: QN = 1−0.15 = 0.85e−. As QN × 24 ≈ 20, this simple calculation agrees with our result, which shows that only 20 out of the 24 N atoms at the surface are protonated. However, only a new FPMD starting from a fully hydroxylated surface plus water should confirm the validity of this estimation. Such a study is currently undergoing and will be published in the near future.

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We acknowledge important discussion with all the members of the SWASSIT group at SBU and BNL, especially with P. B. Allen, J. T. Muckerman, and M. S. Hybertsen. This work is supported by DOE award numbers DE-FG02-09ER16052, DESC0003871, and DE-FG02-08ER46550. This research was carried out (in part) at the Center for Functional Nanomaterials, Brookhaven National Laboratory, which is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, under Contract No. DE-AC02-98CH10886



(1) Maeda, K.; Teramura, K.; Lu, D.; Takata, T.; Saito, N.; Inoue, Y.; Domen, K. Nature 2006, 440, 295−295. (2) Maeda, K.; Teramura, K.; Saito, N.; Saito, N.; Inoue, Y.; Domen, K. Bull. Chem. Soc. Jpn. 2007, 80, 1004−1010. (3) Shen, X.; Allen, P. B.; Hybertsen, M. S.; Muckerman, J. T. J. Phys. Chem. C 2009, 113, 3365−3368. (4) Li, L.; Muckerman, J. T.; Hybertsen, M. S.; Allen, P. B. Phys. Rev. B. 2011, 83, 134202−1−134202−6. (5) Han, W.-Q.; Liu, Z.; Yu, H.-G. Appl. Phys. Lett. 2010, 96, 183112−1−183112−3. (6) Maeda, K.; Domen, K. J. Phys. Chem. Lett. 2010, 1 (18), 2655− 2661. (7) Vittadini, A.; Selloni, A.; Rotzinger, F. P.; Grätzel, M. Phys. Rev. Lett. 1998, 81, 2954−2957. (8) Aschauer, U.; He, Y.; Cheng, H.; Li, S.-C.; Diebold, U.; Selloni, A. J. Phys. Chem. C 2010, 114, 1278−1284. (9) Diebold, U.; Ruzycki, N.; Herman, G. S.; Selloni, A. Catal. Today 2003, 85, 93−100. (10) Shen, X.; Small, Y. A.; Wang, J.; Allen, P. B.; Fernandez-Serra, M. V.; Hybertsen, M. S.; Muckerman, J. T. J. Phys. Chem. C 2010, 114, 13695−13704. (11) Wang, D.; Pierre, A.; Kibria, M. G.; Cui, K.; Han, X.; Bevan, K. H.; Guo, H.; Paradis, S.; Hakima, A.-R.; Mi, Z. Nano Lett. 2011, 11, 2353−2357. (12) Fernandez-Serra, M. V.; Artacho, E. J. Chem. Phys. 2004, 121, 11136−11144. (13) Wang, J.; Roman-Perez, G.; Soler, J. M.; Artacho, E.; FernandezSerra, M.-V. J. Chem. Phys. 2011, 134, 024516−1−024516−10. (14) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (15) Ordejón, P.; Artacho, E.; Soler, J. M. Phys. Rev. B 1996, 53, R10441−R10444. (16) Soler, J. M.; Artacho, E.; Gale, J. D.; A. García, J. J.; Ordejón, P.; Sánchez-Portal, D. J. Phys.: Condens. Matter 2002, 14, 2745−2779. (17) Troullier, N.; Martins, J. L. Phys. Rev. B 1991, 43, 1993,2006. (18) Kleinman, L.; Bylander, D. M. Phys. Rev. Lett. 1982, 48, 1425− 1428. (19) Pamuk, B.; Soler, J. M.; Ramírez, R.; Herrero, C. P.; Stephens, P. W.; Allen, P. B.; Fernández-Serra, M.-V. Phys. Rev. Lett. 2012, 108, 193003−1−193003−5. (20) Fiorentini, V.; Methfessel, M.; Scheffler, M. Phys. Rev. B 1993, 47, 13353−13362. (21) Lilienfeld, O. A. V.; Schultz, P. A. Phys. Rev. B 2008, 77, 115202−1−115202−8. (22) Sulpizi, M.; Sprik, M. Phys. Chem. Chem. Phys. 2008, 10, 5238− 5249. (23) Sprik, M. Chem. Phys. 2000, 258, 139−150. (24) Ivanov, I.; Chen, B.; Raugei, S.; Klein, M. L. J. Phys. Chem. B 2006, 110, 6365−6371. (25) Simon, C.; Ciccotti, G.; Klein, M. L. ChemPhysChem 2007, 8, 2072−2076.



CONCLUSIONS The semiconducting GaN(101̅0) surface has shown to be very reactive, spontaneously dissociating water. The dissociation events happen very fast, and the hydroxilation process at the surface is completed within the first 1.5 ps of simulation. We observed that ∼83% of water molecules at the surface dissociate to OH− and H+, which bind to surface Ga and N, respectively. It was possible to identify two types of hydroxyl ions at the surface, which form dimer pairs with a characteristic structural motif which facilitates the proton diffusion along the surface. This Ga−OH− group is analogous to a water molecule and the Ga−H2O to a hydronium H3O+ with Ga replacing one of the protons. This chemical similarity makes the surface protonation and deprotonation activity of adsorbed surface water molecules very high. In effect, the surface acts as a catalyst for water autodisociation, releasing protons to the bulk liquid region and maintaining the corresponding hydroxyl ions adsorbed. Because of this we were able to compute and effective surface pKa from equilibrium first principles molecular dynamics simulations. The surface is very acidic, with pKa = 3.0 ± 0.1. We have shown that this acidity is in agreement with the observed degree of water dissociation at the surface. We have also shown that the HOMO levels of OH− ions adsorbed at the surface sit above the valence band edge of the semiconductor bulk, allowing them to receive the photogenerated holes. However, to understand which surface species will be the most likely to receive the hole and to initiate the oxidation reaction,10 a full nonadiabatic electron dynamics of the solvated surface is necessary. This is beyond of the scope of this contribution and will be studied in the future, using the structures obtained in this study as our starting point. The facilitated proton diffusion at the surface can be favorable for the transport of the protons from the oxygen evolution reaction site at the semiconductor surface to the hydrogen evolution site in an heterogeneous water splitting photocatalysts. We suggest that this natural acidity might also be one of the reasons behind the enhanced H production activity observed in GaN nanowires with (1010̅ ) terminated surfaces.11



REFERENCES

AUTHOR INFORMATION

Corresponding Author

* E-mail: [email protected]. 14388

dx.doi.org/10.1021/jp302793s | J. Phys. Chem. C 2012, 116, 14382−14389

The Journal of Physical Chemistry C

Article

(26) Sulpizi, M.; Gaigeot, M.-P.; Sprik, M. J. Chem. Theory Comput. 2012, 8, 1037−1047. (27) Leung, K.; Nielsen, I. M. B.; Criscenti, L. J. J. Am. Chem. Soc. 2009, 131, 18358−18365. (28) Cheng, H.; Selloni, A. Langmuir 2010, 26, 11518−11525. (29) Poissier, A.; Ganeshan, S.; Fernandez-Serra, M. V. Phys. Chem. Chem. Phys. 2011, 13, 3375−3384. (30) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397, 601−604. (31) Hassanali, A.; Prakash, M. K.; Eshet, H.; Parrinello, M. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 20410−20415. (32) Berkelbach, T. C.; Lee, H.-S.; Tuckerman, M. E. Phys. Rev. Lett. 2009, 103, 238302−1−238302−4. (33) Agmon, N. Chem. Phys. Lett. 1995, 244, 456−462. (34) Woutersen, S.; Bakker, H. J. Phys. Rev. Lett. 2006, 96, 138305− 1−138305−4. (35) Tielrooij, K. J.; Timmer, R. L. A.; Bakker, H. J.; Bonn, M. Phys. Rev. Lett. 2009, 102, 198303−1−198303−4. (36) Hiemstra, T.; Venema, P.; Riemsdijk, W. J. Colloid Interface Sci. 1996, 184, 680−692. (37) Pashley, M. D. Phys. Rev. B 1989, 40, 10481−10487.

14389

dx.doi.org/10.1021/jp302793s | J. Phys. Chem. C 2012, 116, 14382−14389