Physisorption Structure of Water on the GaN Polar Surface: Force

ARTICLE pubs.acs.org/JPCC

Physisorption Structure of Water on the GaN Polar Surface: Force Field Development and Molecular Dynamics Simulations Osbert Zheng Tan,*,† Michael C. H. Wu,† Viorel Chihaia,‡ and Jer-Lai Kuo*,† † ‡

Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei, Taiwan Institute of Physical Chemistry, Romanian Academy, Bucharest, Romania ABSTRACT: The adsorption mechanism of water on the GaN (0001) polar surface is investigated via both the Density Functional Theory (DFT) method and its derived classical force field. The physisorption binding energy and the adsorption geometry of the water molecule on the clean Ga-terminated surface are analyzed via the first-principle static calculations. The adsorption energy hypersurfaces are then extracted to be used in the fitting of the interaction potentials between water and GaN. Classical molecular dynamics (MD) simulations based on the developed force field are performed for the interfacial system of liquid water and the GaN surface slab. From our computations, the interfacial water exhibits significant oscillatory profiles for the atomic densities and the molecular orientations. Further data analysis suggests a highly confined first layer with the O being locked right upon the surface Ga atoms and the H pointing toward the neighboring O to form the weakened hydrogen bonds. A bilayer configuration with opposite dipole orientations is consequently characterized as the wetting structure on the GaN polar surface and is explained by the anisotropic perturbations from the surface polar sites. Our simulations would be helpful to provide an atomistic picture for the water adsorption configuration on this semiconductor surface and would be useful in the relevant nanofluidic and nanoengineering applications.

I. INTRODUCTION The study of the water/solid interface has been of extensive interest for more than a decade, due to its important applications in wetting, liquid flow, lubrication, phase change, heterogeneous catalysis, etc.1
3 Particularly, for water getting in contact with the solid substances, within a few molecular layers from the interfaces the liquid structure is observed to be significantly rearranged, due to the strong interactions between the water molecules and the wall. Such restructuring of the interfacial water has been evidenced in many interfacial systems, including water on graphene,4 mica,5,6 SiC,7 TiO2,8 Si,9 and NaCl10 surfaces, and they are extensively studied both theoretically and experimentally. Water in the vicinity of solid interfaces is found to have increased viscosity and elasticity due to the surface confinement effect,1 thus exhibiting “solid-like” characteristics that are different from the liquid bulk. For this reason, it has seen many investigations11
15 on the interfacial water properties because of their relevant applications in nanofluidics and nanoengineering. Here, we are focusing on the water/semiconductor interfaces, not only because of the importance of the interfacial water but also because the studies of the water
semiconductor interactions would have their technological implications to understand the semiconductor functionality when the material is exposed to the humid ambient environment. GaN is a material with paramount technological importance, and its broad applications range from optoelectronics and electronics to the gas sensing devices.16,17 The knowledge of the GaN surface characteristics and the surface r 2011 American Chemical Society

interactions with the ambient molecules (particularly water) would therefore be essential for the comprehension of the device operational principles.18,19 Recent first-principle calculations20,21 have identified a dissociative adsorption mechanism for water molecules on the GaN (1010) nonpolar surface, while for the GaN (0001)/(0001) polar surfaces, possible water dissociation mechanisms have been proposed by Chen et al.22 where a moderate reaction barrier is needed for the chemisorption as compared with the spontaneous dissociation of water on the nonpolar surface. The physical wetting structures in the atomistic scale on the GaN surfaces, however, are rarely known, which thus constitutes our current interest. Our theoretical investigations via both DFT and MD simulations are aimed to provide a microscopic model for this interface that can show some interesting water structural profile in the vicinity of the GaN polar surface, and such a model is expected to bridge the gap between the atomic-scale properties of a solid interface and the macroscopic phenomena that can be measured by experiments. Our efforts can also be associated with the recent research on the photocatalytic activities of the GaN (0001) surface23 for the purpose of H2 generation, where the identification of the structural and orientational ordering of the interfacial water molecules would help to understand the mechanism of the photoelectrolysis when water is in contact with the GaN polar surface.24 Received: March 19, 2011 Revised: May 4, 2011 Published: May 23, 2011 11684

dx.doi.org/10.1021/jp202606s | J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C

ARTICLE

Figure 1. (a) Model system for calculations of the adsorption energy hypersurface and the adsorption geometry of water on the GaN bare surface, where the brown represents Ga, the blue N, the white H, and the red O. The passivating H is located at 1 Å below the bottom N. (b) First-principle binding energy hypersurface (symbol data points) of a single water molecule on different adsorption sites of the Ga-terminated surface; the solid lines give the fitted binding energies by using Buckingham potentials for H2O/GaN interactions. dO
Ga indicates the distance between the O and the top Ga plane.

In this paper, we are putting emphasis on the physisorption structure of water on the clean wurtzite GaN (0001) polar surface. By not considering the possible surface chemistry, we are using a nondissociable classical water model interacting with the dipolar GaN Atomic Double Layer (ADL), to explore the thermodynamics and the structured physisorption pattern induced by the strong surface dipolarity. Our current investigations should be considered in parallel with the fact that water on the polar surface is possible to dissociate. Experiment25 suggests that water would chemisorb on the GaN (0001) surface with a sticking coefficient g0.45 and a saturation coverage of 0.46 monolayer, while from the recent transition state calculations,22 only a small reaction barrier of 0.1 eV is required to split water on the Ga-terminated surface. However, the experimentally observed chemisorbed water only covers less than 50% surface sites, and the remaining 50% surface is still exposed to physisorbed water, indicating that the physisorption would be equally important as the chemisorption. In addition, the calculated reaction barrier, even though it is small, can still lead to incomplete dissociations. We are therefore here justifying our motivation by specifically concentrating on the physical structure of liquid water on the defectless GaN surface; the motif is in accordance with the selection of the nondissociable water model and is expected to shed light on the complex surface adsorption mechanisms of this semiconductor. The paper is outlined as below. The second section is dedicated to the computational methods we use for the firstprinciple calculations, force field development, and the molecular dynamics (MD) simulations. It is then followed by the presentation of the results and discussions on the development of a Universal Force Field (UFF) for GaN, the characterization and fitting of water/GaN interactions, and the detailed analysis for the MD on the GaN/water interfacial system. A conclusion is drawn in the end.

II. THEORETICAL METHODS Two model systems are employed for our simulations. For the DFT calculations and the force field parametrizations, a system of a single water molecule on the GaN polar surface slab is adopted, while for the MD simulations, we are using a water slab sandwiched between two GaN polar surfaces with large enough surface areas to probe the possible physisorption patterns at the interfaces.

For the characterization of a single water molecule on the GaN polar surface, we employ Density Functional Theory (DFT) implemented in VASP code.26
29 The projector augmented-wave (PAW) method29 is chosen to represent the pseudopotentials, while the Perdew
Wang 91 functional30 is used for the exchange correlation. For Ga, the 4s, 4p, and 3d electrons are explicitly treated as valence electrons. The electron wave functions are expanded in a plane wave basis set with a cutoff energy of 400 eV. We use the Monkhorst
Pack (MP) scheme31 to generate special k-points, and a (4  4  1) mesh is used for the interface of 1 H2O on the GaN substrate. As shown in Figure 1(a), a system of 6 ADL of GaN in the (0001) orientation is adopted, for which the surface energy is reported to converge at this slab thickness.32 The bottom N face is saturated with pseudo-H atoms where only 0.75e charge is assigned to its valence shell, so that the bulk properties can be converged on the other side of the substrate. The upper layers are terminated with Ga atoms on which water √ molecules can be adsorbed. The surface unit cell is of (2  3) periodicity to achieve the orthorhombicity of the box. The lattice constants are taken from the experiment,33 where a0 = 3.19 Å and c0 = 5.189 Å. The z direction is defined along the surface normal, and the vacuum space is set to be more than 45 Å between two surfaces under the periodic boundary condition. The dipole correction is applied along the z direction to both energy and force,34 so that the artificial electric field in the vacuum space resulting from the dipolarity of the surface slab can vanish.35 To fully capture the interactions within the H2O/GaN interface, we have designed our own force field. The total potential energy of the interfacial system has contributions from the GaN substrate, the water, and the interactions between H2O and GaN Etotal ¼ EGaN þ EH2 O þ EH2 O
GaN

ð1Þ

The potential energy of GaN, EGaN, is designed to follow the functional form of the Universal Force Field (UFF),36 in which a Buckingham potential is used for the van der Waals (VdW) interactions, while the bonded interactions are represented by a harmonic stretching and a nonlinear cosine angular bending, as shown in eq 2. !
 qi qj rij C 1 þ A exp

6 þ kb ðrij
r0 Þ2 EGaN ¼ rij 2 rij F þ kθ ðC0 þ C1 cosðθijk Þ þ C2 cosð2θijk ÞÞ 11685

ð2Þ

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C

ARTICLE

Table 1. Force Field Parameters for the GaN/H2O Interfacial Systema flexible SPC waterb charge (e) Lennard-Jones

O

H


0.82 ε (eV)

0.41 σ (Å)

O
O

0.00674

3.5532

bond stretching

kb (eV Å
2)

r0 (Å)

O
H

48.0595

1

angular bending

kθ (eV rad
2)

θ0 (degree)

H
O
H

3.9695

109.47

The functional is acting on all the four atomic pairs across the interface between O
Ga, O
N, H
Ga, and H
N. The parametrization of the force field for the GaN polar surface and the interactions between H2O and GaN is done via the least-squares fitting procedure implemented in GULP,39 in which an objective function, U, defined as

UFF GaN surface

U ¼

charge (e)

Ga 1.15

N
1.15

Buckingham

A (eV)

F (Å)

C (eV Å6)

Ga
Ga

6260

0.1

0

Ga
N

43.4

0.7

0

N
N

4300

0.12

0

harmonic stretching

kb (eV Å
2)

r0 (Å)

N
Ga

13.2

1.95

nonlinear three-body N
Ga
N

kθ (eV) 10.021

θ0 (degree) 109.47

Ga
N
Ga

10.6755

109.47

GaN
H2O interactions Buckingham

A (eV)

F (Å)

C (eV Å6)

O
Ga

15840.164

0.2

84.852992

N
O

611.33305

0.35

0

H
Ga

9614.862

0.23

0

H
N

1433.2865

0.27

0.003600291

a

By convention, Coulombic interactions are excluded between 1
2 and 1
3 bonded atoms in both water and GaN substrates. b See refs 37 and 38.

The first term gives the Coulombic interaction with qi to be the atomic charge and rij to be the interatomic distance; A, F, and C are parameters in Buckingham potential; kb is the force constant for harmonic bond stretching; r0 is the equilibrium bond length between Ga and N; kθ is the force constant for the nonlinear three-body term; and θijk is the bending angle between two bonds. The parameters C0, C1, and C2 are solely dependent on the equilibrium bending angle θ0, and the functionals (as found in ref 36) are C2 = 1/(4 sin2 θ0), C1 =
4C2 cos θ0, and C0 = C2(2 cos2 θ0 þ 1). It is worth noting that the Lennard-Jones functional presented in the UFF is replaced with the Buckingham here for the VdW interactions. For water, the simple point charge (SPC) model37 with flexible extension38 is implemented. 2 ! !6 3 12 qi qj σ σ 5 þ 1 kb ðrij
r0 Þ2 EH2 O ¼ þ ε4
2 rij rij 2 rij 1 þ kθ ðθijk
θ0 Þ2 2

bending as shown in the last two terms. The flexible SPC parameters are shown in Table 1. The interactions between H2O and GaN are represented by the Coulombic and Buckingham potentials, as shown in eq 4.
 qi qj rij C EH2 O
GaN ¼ þ A exp
ð4Þ
6 rij F rij

ð3Þ

As given in eq 3, the Lennard-Jones functional is used with empirical parameters ε and σ. The intramolecular potential comprises a harmonic bond stretching and a harmonic angular

∑i wi ½fi ðλÞ
Fi 2

ð5Þ

is used to describe the goodness of fit of each individual set. The calculated properties fi are dependent on the parameter set λ. During the fitting, the objective function U is minimized so that fi(λ) can approach their reference values Fi by varying the parameter set λ. The weight of each data item in the fit, wi, determines how well the final fit will reproduce each property. In general, the fitting is done in an iterative manner by starting from some initial guess values, and the fitted parameters are then reinput into the least-squares procedure and repeated and so on until the fitting error is minimized. The GaN UFF is initially parametrized for the wurtzite bulk, based on various bulk properties including lattice constants, bulk modulus, and elastic constants. Both ionic and covalent interactions are involved via the Coulombic and the bonded potential functionals, as compared with the previously developed force field for GaN, where either the ionicity or the covalency of the crystal is considered,40
46 but not for the combination of the two according to our best knowledge. The force field for the GaN polar surface slab is fitted to the surface tetrahedral structure, lattice constants, and the bulk modulus. The parameters of the interaction potential between H2O and GaN are fitted to the geometry and the energy hypersurface obtained by pushing a water molecule toward the GaN (0001) surface where the binding energies as a function of the adsorption distance are collected from the first-principle calculations. More details on the force field development will be given in the following sections. MD simulation is performed in GULP39 by using the parametrized force field, while the electrostatics is calculated via the Ewald summation scheme. As shown in Figure 4(a), the system comprises 680 water molecules with thickness around 42 Å sandwiched between two GaN polar surface slabs, where the number of water molecules is designed to fill in a volume corresponding to the experimental density of liquid water (∼1 g/cm3). Each slab contains 7 ADL along the (0001) orientation, and at both interfaces water molecules are facing the Ga-terminated surfaces. It is worthwhile to mention here that the employment of two surface slabs with opposite net dipoles is to ensure water to be able to fill in a space without the artificial electric field resulted from the dipolarity of the GaN surface slab. We are using an orthorhombic supercell with the surface area designed to be 22.100  22.330 Å2. The bottom one layer of both GaN slabs is fixed to mimic the bulk solid substrate, while the remaining six layers are fully relaxed during the MD and interact with the H2O molecules. The two N-faced pseudosurfaces are separated around 67 Å in the periodic boundary 11686

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C

ARTICLE

Table 2. UFF Parameters for the Wurtzite GaN Bulk, Where the Ga and N Atomic Charges Are Taken from the Kawamura Modela Ga

a

Table 3. Comparison of Calculated and Experimental Results for the GaN Bulk Properties, Where B Represents the Bulk Modulus and Cij Represents the Elastic Constants Matrix

N

bulk properties

calcd

exptla

charge (e)

1.15


1.15

a0 (Å)

3.187

3.19

Buckingham

A (eV)

F (Å)

C (eV Å6)

c0 (Å)

5.207

5.189 0.377

Ga
Ga

6260

0.1

0

u

0.3749

Ga
N

43.4

0.7

0

B (GPa)

236.1

237

N
N harmonic stretching

4300 kb (eV Å
2)

0.12 r0 (Å)

0

C11 (GPa) C12 (GPa)

405.3 138.7

391 143

N
Ga

13.2

1.95

C13 (GPa)

163.1

108

nonlinear three-body

kθ (eV)

θ0 (degree)

C33 (GPa)

384.6

391.2

N
Ga
N

4.021

109.47

C44 (GPa)

126.9

103

Ga
N
Ga

4.6755

109.47

C66 (GPa)

133.3

124

See refs 47 and 48.

condition. The MD is carried out in an NVT ensemble at 300 K. The system is initially relaxed for 0.3 ns using the velocity scaling scheme, and the data are collected in the following 1.7 ns using the Nose-Hoover thermostat. The time step is set to be 0.5 fs.

III. RESULTS AND DISCUSSION A. Development of UFF for GaN and Its Polar Surface. Previous development of the GaN force field has seen Stillinger
Weber40,41 and Tersoff43
45 potentials by focusing on the covalency of the crystal, where the ionic character and the charged polar interactions are completely neglected. Such potentials have been applied to study the crystalline defect and the GaN surface thin film growth where the effect of short-range bonded interactions can be significant. On the other hand, it has seen the development of the partial-charge potential42 and the polarizable core
shell potential46 for GaN by solely focusing on the pairwise ionic interactions. These models can reasonably reproduce the dielectric, elastic properties of the bulk as well as the phase transition mechanisms. In addition, more complex functional form of the two-body potential47,48 (known as Kawamura potential) has been designed for the GaN bulk and is used to investigate the melting point of the wurtzite GaN crystal. For the study of the H2O/GaN interface, the GaN surface polarity plays an important role in prescribing the dipolar interactions with water; thus, the bond order type potentials that ignore the electrostatics would not be suitable in this case. The partial-charge, core
shell, and Kawamura models can well describe the bulk properties as well as the relaxation of the GaN nonpolar surface;46 however, for the polar surface slab, a force field that only considers the two-body ionic interaction would lead to a surface reconstruction from a wurtzite structure to a flat graphitic structure according to our simulations. Such reconstruction, as can be attributed to the destabilizing dipole along the surface normal, however, is artificial, since a Ga-terminated polar surface considered in the experiment is essentially tetrahedral.25,49 We noticed that in first-principle calculations, when similar reconstruction happens, it can simply be healed by increasing the number of layers32 (surface stabilization can be achieved in 6 ADLs), while for the force field, the artifact would persist up to 20 ADLs according to our calculations, indicating that the pairwise potentials are insufficient to create a stable surface structure.

a

See refs 50 and 51.

It is thus appealing to develop a new version of the force field for a GaN polar surface that simultaneously contains two-body and three-body terms to account for the ionicity and the covalency, so that the surface tetrahedrality can be maintained while keeping the electrostatics of the polar films. To do this, we follow the functional form of UFF36 which was used for the tetrahedral Ga and N, except for the VdW term that we have changed to Buckingham potential as introduced in Section II. The new UFF is first parametrized for the wurtzite GaN bulk based on the lattice constants, bulk modulus, and elastic constants with well-referred to experimental values.50,51 As shown in Tables 2 and 3, the calculated lattice constants are in good agreement with the experimental values by having no more than 0.4% deviations. In addition, the bulk modulus and elastic constants are also well reproduced within the uncertainty limits of the experimental reference data.42 The parameters for the GaN polar surface (given in Table 1) are fitted to the tetrahedral surface structure, lattice constants, and the bulk modulus. The parametrized surface force field can give good predictions for the lattice constants (a0 = 3.187 Å, c0 = 5.206 Å, u = 0.375) and the bulk modulus (B = 236.1 GPa) when it is applied on the bulk system, while the elastic constants are observed to have larger deviations from the experimental values. We are making a compromise herein that the force field is primarily designed to maintain the surface tetrahedrality so that the elastic properties are slightly sacrificed. Given this insufficiency, it is worthwhile to mention that, to study the H2O/GaN interface, our current force field would be satisfactory for the purpose of providing a stabilized polar surface that can interact with water. To achieve a better interfacial coupling between the substrate and water, phonon frequencies of the GaN slab would be good surface properties to fit, and it is intended as future work. B. Characterization and Fitting of Water/GaN Interactions. The interaction between water and GaN is characterized by calculating the binding energy of a single water molecule on different adsorption sites of the surface as a function of the adsorption √ distance. Note that here the water is adsorbed on a (2  3) √surface, so it is equivalent to have a water monolayer with 1/2 3 coverage, which is consistent with our aim to study the water-layering structure on the GaN surface. The binding energy is calculated by subtracting the single-water energy and the substrate energy from the energy of the relaxed interfacial system. The water molecule is located on top of the high-symmetry adsorption sites (Ga-top, N-top, and the center of the hexagon) 11687

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C

ARTICLE

Figure 2. (a) Adsorption pathway of H2O on the center adsorption site. (b) Adsorption pathway of H2O on the N-top adsorption site. (c) Adsorption pathway of H2O on the Ga-top adsorption site. Note that the numbers given in the bottom row indicate the distances between O and the Gaterminated plane.

Figure 3. Interatomic potential across the interface. The energy is calculated by using the parametrized force field including both Coulombic and Buckingham potentials for H2O/GaN interactions.

with the vertical-O orientation as the initial geometry (as shown in Figure 1(a) where O is facing the surface and the molecular dipole is pointing along the surface normal). The system is then relaxed by fixing the whole substrate, while for water the O is fixed and only the H atoms are allowed to relax. The adsorption energy hypersurface is obtained by locating the O at different adsorption distances from the surface, while the H atoms are starting from the same initial orientation and optimized to the equilibrium configurations where the binding energies are collected. The optimized water structures show little discrepancy with their respective initial structures for different adsorption sites and different adsorption distances, where the vertical-O orientation is still maintained as the vertical-O one. It can be observed in Figure 1(b) that the Ga-top adsorption site gives the most significant attraction and repulsion to the water molecule, where a binding energy minimum of
0.74 eV is found at 2.01 Å adsorption distance. In contrast, the N-top adsorption site gives a strong repulsion to the water presumably due to the strong electrostatic repulsion between O and N. The adsorption on the center site shows moderate adsorption energy between the Gatop and the N-top ones, while the repulsion is slightly reduced and the water molecule is allowed to approach the surface by

smaller activation energy. For water with vertical-O orientation, all the adsorption energy profiles can be well-categorized into the physisorption class, which actually favors the fitting of the interface interaction with the VdW-type potentials. To investigate how the interfacial interaction affects the water orientation when the molecule is approaching the surface, we perform an adsorption pathway analysis to track the H orientation along the minimal energy pathways of molecules adsorbing on a distinct surface site. As given in Figure 2, the molecule is starting from a vertical-H orientation where one of the H is facing the surface, while the molecular dipole is aligned parallel to the surface (as shown in the rightmost configurations for all adsorption sites in Figure 2). The H is then relaxed by fixing the O at a relatively far distance from the surface, and the structure from the last relaxation is subsequently moved toward the surface by 0.3 Å where the O is fixed again and the H orientations are optimized. The procedure is repeated until the O is as close as 1.87 Å from the top Ga plane. The obtained water geometry on each surface site is therefore continuously optimized along the ground state energy pathway, showing the equilibrium configuration and orientation of a molecule when approaching the interface. Figure 2 displays the optimized structures along the adsorption pathways. All three pathways exhibit that the H is gradually repelled away when the water molecule gets closer to the GaN, while the reason can be tentatively attributed to the steric effect. The development of the force field for the interaction between H2O and GaN is done by fitting the four Buckingham potentials (acting on O
Ga, O
N, H
Ga, and H
N) simultaneously to the three binding energy hypersurfaces obtained with the vertical-O orientation. As given in Figure 1(b), the fitted adsorption energy can well reproduce the first-principle hypersurface on all the adsorption sites. More fitting weight is given to the Ga-top site due to its importance in prescribing the repulsion/attraction interactions with water, while the fitting error (estimated as (fitted data
benchmark data)/(benchmark data)) is kept within 11%. The fitted interatomic potentials across the interface are displayed in Figure 3. C. Molecular Dynamics Simulation for the Water/GaN Interface. The MD simulation is done by putting the water film 11688

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C

ARTICLE

Figure 4. (a) Model system for the MD simulation of water sandwiched between GaN polar surface slabs. Brown, blue, white, and red spheres represent Ga, N, H, and O atoms respectively. (b) Density profile of water along the z axis by using the bin size of 0.0625 Å, where the origin is chosen to be the origin of the unit cell. The error bar is calculated as the standard deviation of the density profiles from ensembles of 100 ps trajectories within the production run; the maximum of the error bar is observed to be 0.08 g/cm3 across the whole z axis.

Table 4. Location of the Peak Positions in Density Profiles, Where The Second Peak in the Total Density Profile Is Estimated As the Shoulder Positiona first peak

second peak

third peak

fourth peak

(Å)

(Å)

(Å)

(Å)

total density O density

2.16 2.16

2.91 2.91

4.22 4.22

5.53 5.53

H density

2.41

4.84

energy

2.01

hypersurface profile on Ga-top a

The last row gives the position of the minimum of the energy hypersurface on Ga-top.

between two Ga-terminated surfaces with the thickness around 42 Å. As shown in Figure 4(b), the density of water is converged to 1 g/cm3 in the middle region with the thickness around 17 Å which mimics the bulk water density, while large fluctuations of the density profile are observed on the two sides due to the interfacial perturbations. The first density peak appearing at 2.16 Å (Table 4) away from the surface reaches the height of 5.47 ( 0.08 g/cm3, demonstrating a very significant rearrangement of the water molecule due to the presence of the surface repulsion/ attraction from GaN. The issue of the density convergence with respect to the bin size is examined, as indicated in Figure 5, where the choices of 0.0625 and 0.125 Å bin widths give rise to indistinguishable density profiles. The magnitude of the first maximum converges into the difference of 0.07 g/cm3 by choosing the above two bin sizes, revealing that our density profile is insensitive to this choice. For discussion purposes, the lower half of the water film on GaN is divided into six regions as shown in Figure 5. The division is based on the total density profile where the boundaries are placed at the local minima of the profile, except for the boundary between regions 1 and 2 which is put at the inflection point (the point where the curvature changes sign) to distinguish the shoulder from the first peak. The total density profile gives four

Figure 5. Density and orientational profile of water on the lower GaN surface, where the origin is chosen to be the mean position of Ga in the top layer during the MD. From top to bottom, it is showing the total water density, O density, H density, and the average molecular dipole orientation as a function of the distance from the Ga-terminated surface. In all the figures, bin size of 0.0625 Å is used, except for the first one where both 0.125 and 0.0625 Å are used to have a convergence test for the density profile. In the bottom figure, φ represents the polar angle between the molecular dipole and the surface normal as illustrated in Figure 6. The boundaries specified by dashed vertical lines between the regions 1
2, 2
3, 3
4, 4
5, and 5
6 are located at 2.78, 3.84, 4.47, 7.09, and 10.41 Å, respectively.

peaks (Table 4) at 2.16, 2.91, 4.22, and 5.53 Å in their corresponding regions; note that the second peak denotes the estimated shoulder position, while the third peak indicates the little hump between the shoulder and the large peak in region 4. The positions of the O density peaks coincide with those in total density profile, whereas the H density profile only gives two remarkable peaks at 2.41 and 4.84 Å within region 1 and 4, respectively. In the first region, the H density peak lying at 0.25 Å above the O density peak suggests a hydrogen-up orientation of the water molecule. The first layer orientation is also suggested by the order parameter profile given in Figure 5, where the average cosine of the molecular dipole polar angle is extracted along the z axis to account for the orientational ordering. As shown in the region 1, almost gives positive value with a maximum of 0.66, demonstrating that the dipole in the first layer is pointing away from the surface. The observation combined with the image that the H density peak locates above the O one is in good agreement with the molecular configurations at around 2 Å from the surface in the adsorption pathway analysis (Figure 2). Nevertheless, there is 0.15 Å discrepancy observed between the first O density maximum and the minimum of the Ga-top energy hypersurface. The reason can be attributed to the orientational effect on the hypersurface, where a more horizontal dipole orientation with larger deviation from the surface normal may bring the equilibrium adsorption position slightly further from the surface than the one determined by the vertical dipole orientation. 11689

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C

Figure 6. Distribution of molecular dipole orientations within each region specified in Figure 5. φ is defined as the angle between the dipole vector of the water molecule and the surface normal.

In region 2, a shoulder can be seen in the total density profile, while a small peak is correspondingly observed in the O density profile. The molecular dipole is generally pointing down toward the surface according to the negative value of in this region (Figure 5). The H density distribution in region 2 is relatively diffuse indicating that there is no in-plane concentration of the hydrogen atoms. Region 3 visualizes a very interesting small hump in the total density and the O density profiles, and is of positive value but with a smaller maximum compared to that in region 1, indicating a dipole upward configuration. In region 4, the second large peak in the total density and the O density profiles is seen, whereas in the H density profile there is a small density gathering in the beginning of the region, which can come from the water molecules with O located at both region 3 and 4. The dipole orientation in this region only sees damped oscillations. In region 5, residue fluctuations are observed in all the density and orientational profiles, while in region 6 the fluctuations converge to their bulk values. Figure 6 shows the molecular dipole orientation distributions within each region. In region 1, the main peak at cos φ = 0.75 corresponds to a dipole orientation directed away from the surface with a preferential polar angle φ = 41.4°, while in region 2, the main distribution of cos φ in the negative value with a peak at cos φ =
0.51 (φ = 120.7°) prescribes a dipole downward configuration. Region 3 sees a sharp distribution at cos φ = 1, demonstrating a large fraction of molecules orientating with their H upward in that small density hump. No significant orientational preference is detected in regions 4 and 5 (the flat distribution in region 4 has the reason that the region contains both dipole-up and dipole-down configurations according to the orientational profile in Figure 5). In region 6, the molecules tend to have equal distributions over all the orientations, which mimics the isotropic liquid bulk. To complement the information from the structural and the orientational profiles along the z axis, we report in Figure 7 the 2D density distributions of O and H atoms in different regions. The O density distribution in region 1 suggests that the oxygen atoms are highly confined within the areas right upon the gallium atoms in the top layer, which is consistent with the first-principle calculations that the water molecule on the Ga-top has the largest binding energy. The first region H distribution shows very interesting patterns that it only covers the areas connecting the

ARTICLE

Ga atoms by not occupying the N-top and the center site. The formation of this particular pattern (as can be visualized more clearly from the average H density distribution in the surface unit cell shown at the bottom of Figure 7) can be attributed to the tendency to form the in-plane hydrogen bond between the horizontally aligned OH and the O located above Ga. The hydrogen bond strength would be relatively weak due to the confinement of O on the Ga-top where the O interatomic distance is found to be around 3.19 Å based on the GaN lattice constant. In the second region, it sees less structural feature in the O density distribution which means more translational freedom of the O in the shoulder region. H atoms however are observed to be highly confined upon Ga atoms. Provided the O is also located in region 2, it can indicate an H-down orientation combined with the fact that the molecular dipole is pointing downward in region 2 (Figure 5). The H pattern in this region can also be attributed to the region 1 molecules with the O confined upon the Ga-top and one of the H vertically protruding into the region 2. In region 3, oxygen atoms adopt some subtle distribution with a hexagonal symmetry (one can refer to the average O density map within the surface unit cell shown at the bottom of Figure 7) around the centers which coincide with the Ga atom positions. The H patterning is not as significant as that in region 2, and it can have the contributions from water with O located at region 3 and both neighboring regions. The delicate hexagonal O pattern and the dipole-up configuration in region 3 probably result from the establishment of a complex hydrogen-bonded network near the surface. From region 4 to region 6, more homogeneous distributions of the O and H atomic density are observed, indicating a reduced surface impact on the anisotropic rearrangement of the molecules. D. Discussion. From the examination of the orientational profile and the spatial density distributions, the first two regions corresponding to the sharp peak and the associated shoulder seem to give very ordered structures. In region 1, the O atoms are located upon the Ga-top with the H slightly tilting away from the surface or even vertically protruding into the next region; while in the shoulder, the H atoms can have the possibility to point downward upon the Ga-top with more translational freedom for the O. The bottom layer arrangement is schematically shown in Figure 8. The dipole up
down configuration seemingly constitutes a dynamic bilayer which can be explained by the anisotropic perturbations from the GaN on the interfacial water structures. It is likely that the molecules first tend to adsorb on the surface polar sites, i.e., Ga-top sites, via the direct polar interactions, so that the molecular dipole is resultantly coaligned with the GaN ADL dipole moment. The molecules beyond the first layer would thus have dipoles with the opposite orientation to compensate the surface dipolarity. To characterize the surface coverage for the molecules in different regions, we also calculate the 2D number density by counting the average number of water molecules (based on the number of O) within each region and dividing the surface area. The surface number density of region 1 is calculated to be 0.083 Å
2, while region 2 is 0.034 Å
2; the full monolayer coverage (1 water adsorbing on 1  1 surface) would correspond to a number density of 0.113 Å
2. Thus, it is interesting to see that the peak region actually does not fully occupy the surface sites, and it only has 0.73 monolayer coverage. The shoulder region is relatively dilute and is of 0.3 monolayer coverage. The combination of the two regions leads to a 1.03 monolayer, indicating that within the film thickness of around 3.84 Å away from the substrate a full surface-bound monolayer 11690

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C

ARTICLE

Figure 7. Normalized O and H density in-plane distribution in different regions. The region index is in accordance with the region division signified in the density profile along the z axis in Figure 5. The density in each region for each atomic species is normalized to the values from 0 to 1 corresponding to the color variations from sharp blue to sharp red, so that the contrast can indicate the spatial concentration of atoms. The mean positions of the top GaN double layer are also plotted with the large white circle representing Ga and the small white circle representing √ N. The figures in the bottom row give the normalized atomic densities of the first three regions averaged onto the GaN surface unit cell with (1  3) periodicity, to highlight the atomic spatial patterning due to the surface impact.

can form, with about 30% of molecules desorbed from the Ga-top sites by having H pointing down to dynamically link to the bottom aligned molecules via the hydrogen bond network. Hereafter, we will still call this structure a virtual dynamic bilayer for the sake of convenience. The surface number density of the third region only reaches 0.014 Å
2, corresponding to a coverage of 0.12 monolayer. It indicates the formation of a very dilute layer at the transition region between the first and the second major density peaks. The slightly gathered O atoms distribute on a subtle hexagonal pattern by not occupying the Ga-top sites, which can be tentatively attributed to the formation of the extended hydrogen bond network with the surface contact bilayer. The dipole-up configuration can be explained by the induced orientational ordering due to the residue surface polarity from the contact bilayer structure, as evidenced by the damped oscillations in the orientational profile (Figure 5). In addition, it is worthwhile to mention that the currently identified water structure on the surface of GaN (including a bilayer and a third layer with induced ordering) is completely based on a nondissociable water model interacting with the polar GaN surface with defectless tetrahedral structure. Such an

unreconstructed clean surface is purported by Sung et al.49 from their low-energy electron diffraction (LEED) studies, but it is in contradiction with many other experimental results52
56 that the GaN (0001) surface exhibits a large variety of reconstructions, such as (1  1), (2  2), (4  4), (5  5), and (6  4). The reconstructed surface is reported18 to be highly reductive thus easy to grab oxygen from the adsorbing molecular species. The interaction of water with the clean unreconstructed surface however remains ambiguous except for one article25 reporting that H2O would adsorb dissociatively on the GaN surface with a sticking coefficient g0.45 and a saturation coverage of 0.46 monolayer; however, the provided surface for the adsorbate is an oxidized surface with the residue oxygen coverage in the range of 0.1 monolayer. It is also interesting to see the recent calculations from Chen et al.22 where the water molecule can dissociate on the clean GaN (0001) surface with a moderate reaction energy barrier, but the relevancy of their computations regarding the surface chemistry at finite temperature concerning the reconstruction issues of the GaN is still difficult to justify according to their static transition state calculations. The proposed dynamic bilayer structure based on our current analysis can therefore be used to complement the understanding of the possible 11691

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C

Figure 8. MD snapshot showing molecular configurations of the bottom layers within around 4.82 Å from the surface. (a) Top view. (b) Side view. The number indicates the region divisions according to the specification in Figure 5. In both graphs, the blue dashed lines represent the hydrogen bonds.

chemisorption mechanism, provided relevant experiments can be implemented to elucidate the physical wetting structure of water on a clean unreconstructed surface. It also should be noted that if the chemisorption is taken into account the proposed water structuring at the interface can be disturbed due to the hydrogen bond interactions with the surface hydroxyl and the chemisorbed oxygen.

IV. CONCLUSIONS We have designed a modified UFF for the GaN polar surface that can provide stabilized tetrahedral surface structure with full electrostatics. The force field for the interactions between water and GaN is implemented with Buckingham potentials where the parameters are fitted to the binding energy hypersurface of water adsorbed on different surface sites. MD simulations are performed using the parametrized force field for the water/GaN interface; the results give rise to a highly structured profile for water from the immediate vicinity of the surface up to a few molecular dimensions. The anisotropic perturbations to water from the GaN polar surface sites can be characterized via the sharp water density peak at 2.16 Å associated with a shoulder at 2.91 Å. The molecular dipole in the peak region prefers to orient away from the surface with the preferential tilt angle φ = 41.4°, while in the shoulder region the dipole tends to orient toward the surface with the preferential tilt angle φ = 120.7°. The 2D density distribution of O in the peak region shows a highly confined pattern right upon the Ga-top sites, with H only covering the areas connecting the Ga atoms due to the formation of the weak hydrogen bonds, while H in the shoulder region displays a confined distribution upon the Ga-top which mimics that of O

ARTICLE

in the first region. O in this region however shows a relatively diffuse pattern with more translational freedom. The observed structuring in the contact bilayer agrees well with the DFT calculations in terms of the adsorption distance, adsorption geometry, and adsorption pattern, while the opposite dipole orientations within the bilayer can be easily understood via the strong polar interactions of molecules from region 1 with the surface sites and the dipolar compensation of molecules in region 2. The detected small hump in the density transition region indicates the slight gathering of the O atoms, which however also gives rise to a hexagonal O distribution pattern. The reason is attributed to the induced ordering from the residual polarity and the extended hydrogen bond network with the contact bilayer structure. In summary, we have visualized a highly ordered profile for water on the GaN polar surface. The interfacial structuring of molecules is comparable with that of water on other hydrophilic surfaces (i.e., silicon9 and hydroxylated silica15) in terms of the density and the orientational profile oscillations. The surface dipolarity has led to a significant rearrangement for the interfacial molecules, and the impact is reduced at around 4.47 Å from the surface beyond which the molecular orientational ordering and the spatial patterning are found to be less remarkable. Other dynamic properties, such as the hydrogen bond lifetime and the reorientational relaxation, which can distinguish the interfacial water from the bulk, would be discussed in detail in the forthcoming papers. Our results presented here would be helpful to understand the physical wetting structure of water on this semiconductor polar surface and also provide a basis to investigate the molecular-scale phenomena relevant to nanofluidics and other nanotechnological applications.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]; [email protected]. sinica.edu.tw.

’ ACKNOWLEDGMENT This work is supported in part under Academia Sinica Research Program on NanoScience and Nano Technology and National Science Council (NSC98-2113-M-001-029-MY3) of Taiwan. ’ REFERENCES (1) Antognozzi, M.; Humphris, A. D. L.; Miles, M. J. Appl. Phys. Lett. 2001, 78 (3), 300. (2) Fenter, P.; Sturchio, N. C. Prog. Surf. Sci. 2004, 77 (5
8), 171. (3) Hussain, M.; Anwar, J. J. Am. Chem. Soc. 1999, 121 (37), 8583. (4) Cicero, G.; Grossman, J. C.; Schwegler, E.; Gygi, F.; Galli, G. J. Am. Chem. Soc. 2008, 130 (6), 1871. (5) Fukuma, T.; Ueda, Y.; Yoshioka, S.; Asakawa, H. Phys. Rev. Lett. 2010, 104 (1), 016101. (6) Leng, Y. S.; Cummings, P. T. Phys. Rev. Lett. 2005, 94 (2), 026101. (7) Cicero, G.; Grossman, J. C.; Catellani, A.; Galli, G. J. Am. Chem. Soc. 2005, 127 (18), 6830. (8) Koparde, V. N.; Cummings, P. T. J. Phys. Chem. C 2007, 111 (19), 6920. (9) Xu, D. Y.; Leng, Y. S.; Chen, Y. F.; Li, D. Y. Appl. Phys. Lett. 2009, 94 (20), 201901. 11692

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693

The Journal of Physical Chemistry C (10) Liu, L. M.; Krack, M.; Michaelides, A. J. Chem. Phys. 2009, 130 (23), 234702. (11) Gordillo, M. C.; Marti, J. Phys. Rev. B 2008, 78 (7), 075432. (12) Marti, J.; Gordillo, M. C. Phys. Rev. B 2001, 63 (16), 165430. (13) Marti, J.; Nagy, G.; Gordillo, M. C.; Guardia, E. J. Chem. Phys. 2006, 124 (9), 094703. (14) Marti, J.; Nagy, G.; Guardia, E.; Gordillo, M. C. J. Phys. Chem. B 2006, 110 (47), 23987. (15) Argyris, D.; Tummala, N. R.; Striolo, A.; Cole, D. R. J. Phys. Chem. C 2008, 112 (35), 13587. (16) Pearton, S. J.; Ren, F.; Wang, Y. L.; Chu, B. H.; Chen, K. H.; Chang, C. Y.; Lim, W.; Lin, J. S.; Norton, D. P. Prog. Mater. Sci. 2010, 55 (1), 1. (17) Lee, D. S.; Lee, J. H.; Lee, Y. H.; Lee, D. D. Sens. Actuator, B: Chem. 2003, 89 (3), 305. (18) Lorenz, P.; Gutt, R.; Haensel, T.; Himmerlich, M.; Schaefer, J. A.; Krischok, S. Phys. Status Solidi C 2010, 7 (2), 169. (19) Stutzmann, M.; Steinhoff, G.; Eickhoff, M.; Ambacher, O.; Nebel, C. E.; Schalwig, J.; Neuberger, R.; Muller, G. Diamond Relat. Mater. 2002, 11 (3
6), 886. (20) Shen, X.; Allen, P. B.; Hybertsen, M. S.; Muckerman, J. T. J. Phys. Chem. C 2009, 113 (9), 3365. (21) Shen, X. A.; Small, Y. A.; Wang, J.; Allen, P. B.; Fernandez-Serra, M. V.; Hybertsen, M. S.; Muckerman, J. T. J. Phys. Chem. C 2010, 114 (32), 13695. (22) Chen, P. T.; Sun, C. L.; Hayashi, M. J. Phys. Chem. C 2010, 114 (42), 18228. (23) Fujii, K.; Karasawa, T. K.; Ohkawa, K. Jpn. J. Appl. Phys., Part 2 Lett. Express Lett. 2005, 44 (16
19), L543. (24) Kanemaru, H.; Nosaka, Y.; Hirako, A.; Ohkawa, K.; Kobayashi, T.; Tokunaga, E. Opt. Rev. 2010, 17 (3), 352. (25) Bermudez, V. M.; Long, J. P. Surf. Sci. 2000, 450 (1
2), 98. (26) Kresse, G.; Furthmuller, J. Phys. Rev. B 1996, 54 (16), 11169. (27) Kresse, G.; Furthm€uller, J. Comput. Mater. Sci. 1996, 6 (1), 15. (28) Kresse, G.; Hafner, J. Phys. Rev. B 1993, 47 (1), 558. (29) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59 (3), 1758. (30) Perdew, J. P.; Wang, Y. Phys. Rev. B 1992, 45 (23), 13244. (31) Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13 (12), 5188. (32) Freeman, C. L.; Claeyssens, F.; Allan, N. L.; Harding, J. H. Phys. Rev. Lett. 2006, 96 (6), 066102. (33) Xu, Y.-N.; Ching, W. Y. Phys. Rev. B 1993, 48 (7), 4335. (34) Neugebauer, J.; Scheffler, M. Phys. Rev. B 1992, 46 (24), 16067. (35) Krukowski, S.; Kempisty, P.; Strak, P. J. Appl. Phys. 2009, 105 (11), 113701. (36) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114 (25), 10024. (37) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction models for water in relation to protein hydration. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel: Amsterdam, 1981; p 331. (38) Teleman, O.; J€onsson, B.; Engstr€om, S. Mol. Phys. 1987, 60 (1), 193. (39) Gale, J. D. J. Chem. Soc., Faraday Trans 1997, 93 (4), 629. (40) Aïchoune, N.; Potin, V.; Ruterana, P.; Hairie, A.; Nouet, G.; Paumier, E. Comput. Mater. Sci. 2000, 17 (2
4), 380. (41) Bere, A.; Serra, A. Phys. Rev. B 2002, 65 (20), 205323. (42) Gao, F.; Devanathan, R.; Oda, T.; Weber, W. J. Nucl. Instrum. Methods Phys. Res., Sect. B 2006, 250, 50. (43) Moon, W. H.; Hwang, H. J. Phys. Lett. A 2003, 315 (3
4), 319. (44) Nord, J.; Albe, K.; Erhart, P.; Nordlund, K. J. Phys.: Condens. Matter 2003, 15 (32), 5649. (45) Nord, J.; Nordlund, K.; Keinonen, J.; Albe, K. Nucl. Instrum. Methods Phys. Res., Sect. B 2003, 202, 93. (46) Zapol, P.; Pandey, R.; Gale, J. D. J. Phys.: Condens. Matter 1997, 9 (44), 9517. (47) Harafuji, K.; Tsuchiya, T.; Kawamura, K. Phys. Status Solidi C 2003, 0 (7), 2420. (48) Harafuji, K.; Tsuchiya, T.; Kawamura, K. J. Appl. Phys. 2004, 96 (5), 2501.

ARTICLE

(49) Sung, M. M.; Ahn, J.; Bykov, V.; Rabalais, J. W.; Koleske, D. D.; Wickenden, A. E. Phys. Rev. B 1996, 54 (20), 14652. (50) Polian, A.; Grimsditch, M.; Grzegory, I. J. Appl. Phys. 1996, 79 (6), 3343. (51) Xia, H.; Xia, Q.; Ruoff, A. L. Phys. Rev. B 1993, 47 (19), 12925. (52) Smith, A. R.; Feenstra, R. M.; Greve, D. W.; Shin, M. S.; Skowronski, M.; Neugebauer, J.; Northrup, J. E. J. Vac. Sci. Technol. B 1998, 16 (4), 2242. (53) Smith, A. R.; Feenstra, R. M.; Greve, D. W.; Shin, M. S.; Skowronski, M.; Neugebauer, J.; Northrup, J. E. Surf. Sci. 1999, 423 (1), 70. (54) Timon, V.; Brand, S.; Clark, S. J.; Gibson, M. C.; Abram, R. A. Phys. Rev. B 2005, 72 (3), 035327. (55) Widmann, F.; Daudin, B.; Feuillet, G.; Pelekanos, N.; Rouviere, J. L. Appl. Phys. Lett. 1998, 73 (18), 2642. (56) Yu, Z. X.; Tong, S. Y.; Xu, S. H.; Ma, S.; Wu, H. S. Surf. Rev. Lett. 2003, 10 (6), 831.

11693

dx.doi.org/10.1021/jp202606s |J. Phys. Chem. C 2011, 115, 11684–11693