Molecular surface structure of ice(0001): dynamical low-energy electron diffraction, total-energy calculations and molecular dynamics simulations

surface science ELSEVIER

Surface Science 381 (1997) 190 210

Molecular surface structure of ice(0001 )" dynamical low-energy electron diffraction, total-energy calculations and molecular dynamics simulations N. Materer a,1, U. Starke a,2, A. Barbieri a,3, M.A. Van Hove a,., G.A. Somorjai a, G.-J. Kroes b C. Minot c a Materials Sciences Division, Lawrence Berkeley National Laboratory, University of Cal([ornia, and Department o f Chemistry, University o f CaliJbrnia, Berkeley, CA 94720, USA b Theoretical Chemistry, Department o['Chemistry, Vrije Universiteit, De Boelelaan 1083, 1081 H V Amsterdam, Netherlands Laboratoire de Chimie Thdorique, URA 560 Tour23-22 114, 4 Place Jussieu, Universiff~Pierre et Marie Curie, F-75230 Paris, France

Received 19 March 1996; accepted for publication 9 January 1997

Abstract A structural study of the surface of an ultrathin ice film grown on a Pt( l I l ) substrate was performed using dynamical low-energy electron diffraction (LEED) at 90 K, together with total-energy calculations and molecular dynamics (MD) simulations. This ice film exhibits the common hexagonal phase Ih and exposes the (0001) surface without reconstruction. The surface is terminated as a full-bilayer that maximizes the number of surface hydrogen bonds as confirmed by our total-energy calculations. Both LEED and MD simulations find that the outermost water molecules have enhanced vibrational amplitudes making them practically undetectable by LEED even at 90 K. MD simulations of the half-bilayer terminated surface yield results inconsistent with the LEED findings, thus excluding this model. © 1997 Elsevier Science B.V. Keywords: Ab initio quantum chemical methods and calculations; Auger electron spectroscopy; Crystallization; Electron solid

diffraction; Electron solid interactions, scattering, diffraction; Epitaxy; Low energy electron diffraction (LEED); Low index single crystal surfaces; Molecular dynamics; Single crystal surfaces; Surface melting; Surface relaxation and reconstruction; Surface structure, morphology, roughness, and topography; Vibrations of adsorbed molecules; Water

1. Introduction Ice p l a y s a p e r v a s i v e r o l e in m a n y d a i l y p h e n o m e n a , p a r t i c u l a r l y t h r o u g h its surface. F o r e x a m ple, the u n u s u a l slipperiness o f ice is p o s s i b l y d u e

to s u r f a c e p r e m e l t i n g o r m e l t i n g o f the s u r f a c e b e l o w the b u l k m e l t i n g p o i n t (for reviews, see Refs. [1,2]). Ice is t h o u g h t to c a t a l y z e c h e m i c a l r e a c t i o n s in the E a r t h ' s o z o n e l a y e r a n d it c a n p l a y a role in c l o u d - s e e d i n g a n d t h u n d e r s t o r m s .

* Corresponding author. Fax: + 1 510 4864995; e-mail: [email protected] i Present address: JILA, University of Colorado, Boulder, CO 80309, USA. 2 Present address: Lehrstuhl for FestkOrperphysik Universit~t Erlangen-N~irnberg, Staudtstr. 7, D-91058 Erlangen, Germany. 3 Present address: FEA Inc., 2484 Shattuck Ave., Suite 225, Berkeley, CA 94704, USA. 0039-6028/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. P H S0039-6028 (97)00090-3

N. Materer et al. / Surface Science 381 (1997) 190~10

Previous experimental and theoretical investigations of the surface of bulk ice have focused on adhesion and friction [3], on surface premelting [ 1,2,4-9], on vaporization mechanisms [ 10] and on the adsorption of non-water molecules on ice clusters and amorphous ice [11,12]. In other studies, measurements of O - H vibrations in ice clusters and amorphous ice can identify hydrogen atoms not participating in hydrogen bonds and thus characterize the surface [13]. Most studies so far were performed at temperatures well above 240 K and report the presence of a liquid or quasiliquid layer on ice. Those studies that went below this temperature do not suggest a liquid-like layer. In one particular study, X-ray crystal truncation rod scattering was used to determine the possibly preparation-dependent surface morphology of the (0001) surface of the common ice phase Ih at 266 K, but the surface crystallography [6] was not determined. In this paper, we examine the ice surface both experimentally, using low energy electron diffraction (LEED), and theoretically, using both molecular dynamics and total energy computations, in order to provide a structural model for the surface. One question which is also addressed in this paper is to what extent ice surfaces may show either disordering behavior or larger amplitudes of vibrations at temperatures as low as 90 K. At this temperature, surface melting is not expected to occur. Most experiments investigating surface melting of ice put the temperature at which a liquid-like layer may exist on the ice at no lower than about 243 K ( - 3 0 ° C ) . At low temperatures, one still might expect larger amplitudes of motion and a greater extent of disordering in the surface layer, due to disruption of the hydrogen bonding network at the surface. In the bulk, an H 2 0 molecule is held in place by four hydrogen bonds, whereas, a surface molecule, in its ideal bulk-like arrangement, can form no more than three hydrogen bonds. The issues of surface disorder and enhanced vibrational amplitudes are especially important for structural determinations of molecular surfaces where the bonding is much weaker than in metals or semiconductors. U H V studies of water adsorption on single

191

crystal surfaces have utilized a variety of surface science techniques [ 14]. Electron stimulated desorption ion-angular distribution (ESDIAD) and vibrational spectroscopy have provided orientational information for submonolayer to saturation coverages of water molecules. Experiments on Pt(111 ) using high-resolution electron energy loss spectroscopy ( H R E E L S ) and more recently with infrared reflection absorption spectroscopy (IRAS) detect clustering of the water molecules even in the submonolayer regime [15]. On Ru(0001), a detailed structural study of water adsorption has been conducted [16]. This study found that deuterated water molecules form a single ordered bilayer, which fits the substrate 2D lattice in a (~¢/3x~f3)R30 ° pattern and that the oxygen positions could be determined with automated tensor LEED. In their final structure, the outermost water molecules, which are hydrogenbonded to the molecules underneath, were found to be almost coplanar with the innermost metalbonded molecules. The earliest MD simulations investigating surface disordering of ice studied the melting of small crystallites in the temperature range between 150 and 300 K [17]. Melting was observed to start at the surface. It was argued that the disruption of the hydrogen bonding network that occurs at the surface sets off the surface melting process. Subsequently, the melting of a slab of ice was investigated in somewhat more detail for the temperature range between 190 and 250 K [5]. The more recent simulations show evidence for partial disordering at a temperature as low as 190 K. MD simulations are presently used also to study ice growth [18]. Very recently, the M D technique has also been applied to the simulation of a platinum-water interface [ 19, 20]. In this paper, we will first briefly describe our thin film approach, used to study the bulk ice surface, in Section 2. Provided in Section 3 is a short presentation of our theoretical methods, while Section 4 contains the LEED results. In this section a detailed description of the various surface models and their relative agreement with the experimental intensity vs voltage curves ( L V curves) can be found. Because it was not possible to

192

N. Materer et al. / Sur/ace Science 381 (1997) 190-210

distinguish between models which result in an apparent half-bilayer vs a full-bilayer termination, we address these issues with both total energy computations and molecular dynamics simulations. Section 5 describes total energy computations used to obtain the relative energy differences between various surface models. To investigate the possible influence on the LEED experiments of either a disordered top layer (due to thermal motions which result in the loss of lateral order) or enhanced vibrational amplitudes of the surface molecules, molecular dynamics (MD) simulations of the basal plane or (0001) face of ice at the experimental temperature ( T = 90 K ) were carried out. We simulate both the full- and the half-bilayer termination. The full-bilayer termination model with enhanced vibrational amplitudes is supported by both the total energy computations and MD simulations. Section 7 refines the LEED analysis of this model. Finally, our conclusions are summarized in Section 8.

2. The LEED experiment The LEED experiments were carried out in a stainless steel ultrahigh vacuum chamber (described in Ref. [21]) which was kept at or below 10 10 Torr throughout the experiments. An exposure of 20 L of water vapor at 140 K resulted in the ordered six-fold symmetric LEED pattern shown in Fig. 1 [21,22]. On the P t ( l l l ) surface, no ordered water adsorption structures are formed with the exception of a multilayer ice film [21,22]. Comparing both the LEED spots and the background intensity of the LEED pattern shown in Fig. 1 to the LEED pattern obtained from the clean P t ( l l l ) surface (not shown), one observes both the relatively larger spot size and higher background intensity in the ice film LEED pattern. These features indicate a larger degree of disorder than that found for the clean metal surface. In particular, the background could be indicative of diffuse thermal scattering due to surface vibrations. The LEED patterns measured from the ice film showed no contributions from the P t ( l l l ) substrate; also, the 2D lattice constant computed from

i'i

Fig. 1. The L E E D pattern of the ice film grown on the Pt( I 11 ) surface at 48 eV and normal incidence.

the LEED patterns matches the lattice constant of bulk ice within the LEED uncertainty of a few percent. A lower limit estimate of the film thickness can be obtained from the electron attenuation length 2, which, in turn, can be estimated by the examination of the peak widths in the LEED L V curves. From these considerations we have determined that the thickness of the film, as seen by the LEED electron, must be at least 10 A [22]. After the preparation of the ice film, the crystal was cooled to 90 K. A digital LEED detector [23] was used to measure the impact position of single diffracted electrons and digitally generate the diffraction pattern. We exclude the influence of electron-beam-induced damage on our diffraction patterns, a definite possibility with water molecules and ice [14], since our experiment used an extremely weak picoampere electron beam current that deposits only about 10-2 electrons per experiment per surface unit cell [23]. Intensity vs energy ( ~ V ) curves were measured for 18 beams from 30 to 150 eV. The symmetrically-equivalent beams in three rings of spots in the six-fold symmetric pattern were averaged, resulting in three nonequivalent beams, the (1,0), (1,1) and (2,0), for use in the LEED analysis. The collection of higher-

N. Materer et al. / Surface Science 381 (1997) 190 210

energy I - V data was prevented by the rapid fall of intensity with increasing energy. This rapid fall of intensity was also reproduced in the theoretical I - V curves and is due to the decreasing crosssection for backscattered electrons with increasing energy. The curves were repeatedly measured after four separate surface preparations and found to be reproducible.

3. Theory 3.1. Low-energy electron diffraction theory The LEED simulations utilized the well-established automated tensor LEED method [24]. We modeled the ice surface as a semi-infinite ice sample using composite bulk and surface bilayers which were combined using renormalized forward scattering [24]. For ice, seven scattering phase shifts were computed utilizing a potential obtained, within the muffin-tin approximation, from overlapping atomic wave functions on an Ih ice lattice. Inelastic damping of the electrons was represented by the imaginary part of the inner potential of - 4 eV. Thermal vibrations were represented by temperature dependent phase shifts based on a Debye-Waller model. The rms vibrational amplitudes used to construct the temperature dependent phase shifts were fit to experiment and will be discussed along with the LEED results.

3.2. Total energy theory Total energy calculations were performed at the Hartree Fock level utilizing CRYSTAL [25]. This program carries out total energy calculations for finite slabs made up of two-dimensional layers which are of infinite extent. We have used two basis sets: the 631G* and the 631G** [26 29]. The main difference between these two basis sets is that the G31G** allows for some polarization of each H atom's charge density by the addition of a p-orbital function. All energies reported in this paper are energy gains with respect to the free molecule(s). For the free water molecule, these basis sets give OH distances of 0.947 and 0.943 A,

193

while for the HOH angles they give 105.5 ° and 106 °. These values are to be compared to the gas values of 0.9572 A and 104.52 °. The energy for an intermolecular bond, 5.5 to 5.6kcalmo1-1, is larger than the SCF limit [30,31] (the limit obtained with an infinite basis set), by 1.9 kcal mol 1. Most of the published total energy calculations concern the proton-ordered phases that are stable at high pressure and low temperature [32,33]. The influence of the oxygen coordinates and of the hydrogen positions in the lattice on the total energy are not known; thus, we perform geometry optimizations on several different slab geometries to estimate these effects.

3.3. Molecular dynamics theory Ice surface simulations were carried out in essentially the same manner used in a previous study [5] for a higher range of temperatures (between 190 and 250 K). In the present work, we simulate an ice surface at 90 K without modeling the interaction with the platinum underneath. This seems appropriate as the ice layer we attempt to simulate is thicker than 10 A. This means that the ice layer may be thick enough not to be affected any more by the ice-platinum interface. In these MD simulations, the (0001) surface is modeled using 12 layers of "moving" H20 molecules which form a slab of thickness 21.83 A. The moving molecules are held in place by four layers of fixed H20 molecules underneath. An infinite surface is simulated by replicating the " M D box" containing the H20 molecules in the directions parallel to the surface, using periodic boundary conditions. The MD box has the dimensions x=22.35 A and y=23.23 A, parallel to the surface, and contains 480 H20 molecules, 360 of which move according to Newton's equations of motion. The simulations are started from a configuration generated using a scheme suggested by the work of Cota and Hoover [34]; this configuration has a zero dipole moment and obeys the "ice rules" [35]. We use the TIP4P potential model [36,37] for the interaction between the H20 molecules. An advantage of this potential is that it yields stable hexagonal ice [38]. In the TIP4P model, the water

194

N. Materer el al. / Sur[ace Science 381 (1997) 190 210

molecule is treated as a rigid rotor, disallowing intramolecular vibrational motion. Additional discussions of the TIP4P potential applied to the surface of ice can be found in Ref. [5]. As was previously done [5], the potential is switched off at long range (between 9.5 and 10A for the switching function); this is because, for a slab of molecules in a semi-infinite lattice, it would be both difficult and computationally expensive to calculate the long range interaction using an Ewald-like scheme [39,40]. Furthermore, MD simulations [18] of a slab of ice which did use an Ewald-like scheme [41,42] to evaluate the electrostatic forces, but employed a smaller unit cell than our simulations, yield similar results for the amplitudes associated with the surface phonons. In the present work, these amplitudes are the most important result from the MD simulations. In the simulations, the temperature (90 K) is imposed using a computational analog of a thermostat [43] during an initial 10 ps run. The equations of motion are solved using an improved leapfrog algorithm {44], with a timestep of 1 fs. After imposing the temperature of 90 K, the surface was propagated in time in four runs all lasting 20 ps. The results presented below are (converged) averages of the last run. Comparisons with results using 16 (instead of 12) layers of moving water molecules showed that the results presented below are also converged with respect to the number of layers used in the simulations. The rms amplitude of motion of the water molecules perpendicular to the surface was calculated by computing, for each layer, the average value of

,~(z) = V ( Z o ) - (Zo)'~,

(1)

where Z o is the instantaneous Z-coordinate associated with the oxygen atom of a water molecule and Z is the coordinate perpendicular to the surface. The bar represents an average over the simulation(s). To assess the disorder associated with the vertical displacements of the molecules, density profiles were made by partitioning the sample along the direction perpendicular to the surface into bins 0.1 A wide, and summing the

number of oxygen atoms in each bin at each time step. To monitor disorder associated with translational motion parallel to the surface, a translational order parameter (ST) was calculated by averaging the three order parameters given by Eq. (2) in Ref. [5]. This order parameter measures to what extent the O-atoms in a given layer remain in their ideal lattice positions parallel to the surface (relative to one another). For a perfectly ordered surface, S T = l , while for a perfectly disordered (liquid-like or amorphous) surface, ST is equal to one divided by the number of molecules per layer (30 in the present calculation). Rotational disorder was monitored quantitatively by calculating a rotational order parameter (S,) (Eq. (3) of Ref. [5]). The order parameter SR measures to what extent the orientations of the dipole vectors of the H20 molecules in a given layer deviate from the ideal orientations, of which there are six for each layer. For a perfectly ordered solid, SR = 1, while for a liquid SR =0. More qualitatively, rotational disorder was calculated by calculating distributions ("densities") of cos 0 and ~b, where 0 and ~b are the polar and azimuthal angles that the dipole vector of a water molecule makes with respect to the surface normal. The distributions were calculated by partitioning cos 0 into 100 bins and ~b into 200 bins, and by summing the number of cos 0 and ~b of each water molecule in each bin at each timestep. The simulations of the (0001) surface with the half-bilayer termination were done in essentially the same manner, except that now the upper halfbilayer of molecules was removed before performing the simulations. A slightly longer equilibration time was used; the results presented here are the averages of the last run of six consecutive runs which lasted 20 ps each.

4. Preliminary LEED analysis 4.1. Surface models to be tested

The LEED pattern only fits the two simplest phases among the many different bulk ice phases [45,46]: the unreconstructed (0001)-( 1 × 1 ) surface

N. Materer et al. / Surlhce Science 381 (1997) 190 -210

of hexagonal ice Ih, the most commonly occurring form, and the ( 111 )-( 1 x 1 ) surface of cubic ice Ic, a metastable phase that has been obtained by condensation on metals at about 140 K [14]. Both bulk ice Ih and bulk ice Ic consist of bilayers of tetrahedrally bonded HzO molecules connected by hydrogen bonds. The O layers have hexagonal symmetry and the O planes are spaced by about 2.76 A between bilayers and 0.92 A within bilayers. The O atoms in Ic form the diamond lattice while those in Ih form the corresponding hexagonal lattice shown in Fig. 2. The O layers have an ... ABBAABBA ... stacking sequence in ice Ih, and an ... ABBCCAAB ... sequence in ice Ic, so the distinction only appears after three or four O layers, or about 7 A under the ideally-terminated surface. Each oxygen is surrounded by four other hydrogen-bonded oxygens about 2.76 A away at RT and 2.74 A at 90 K. The H atoms are located between adjacent O atoms in the hydrogen bond. but are off-set asymmetrically by about 0.38 from the midpoint in a largely random fashion and are most likely slightly off-axis [45]. Each O atom is covalently bonded to two nearby H atoms and weakly bonded to two other, more distant H atoms. We have investigated models based on two different terminations of the ice lh(0001) and I c ( l l l ) surfaces, as shown in Fig. 2 for the Ih(0001) model. In the full-bilayer termination model each outermost O atom is hydrogen-bonded to 3 neighboring O atoms, while in the half-bilayer termination model each is hydrogen-bonded to

Full-Bilayer Termination

tlalf-Bilayer Termination

V

Layer 1 Layer 2 Layer 3 2.74/~

""

'

;

Fig. 2. Perspective grazing view of the ice Ih(0001 ) surface with two ideal terminations: full-bilayer termination at left, and halfbilayer termination at right. Large circles are oxygens, with some hydrogens included as small circles with assumed bulklike positions and r a n d o m n e s s to form HzO molecules (covalent bonds are drawn thick, hydrogen bonds thin).

195

only one neighboring O atom. For the full-bilayer lh and Ic models, we also allowed the outermost O to have an enhanced vibrational amplitude. As discussed below, we found that a large enhancement is needed, such that the Debye-Waller factor in back scattering makes the outermost O atoms practically "invisible" (there remains undiminished narrow forward scattering which has a small but detectable effect through multiple scattering paths). Enhanced vibrational amplitudes were not tried for the half-bilayer model because this feature, due to the Debye Waller factor, would effectively, for large amplitudes, make the model equivalent to the full-bilayer termination model. In Section 7.2, we will discuss the optimization of these enhanced vibrational amplitudes. In an attempt to simulate a surface with a mixture of full- and half-bilayer termination, an incoherent average of the scattered intensities resulting from the two terminations was also investigated. In addition, a collapsed bilayer was tried, in which the water molecules forming the outermost bilayer are nearly coplanar. For all models, domain averaging was applied to recover the six-fold symmetry observed in the LEED pattern. For each model, our calculations explored different lateral lattice constants from 4.33 to 4.99 A parallel to the surface the bulk ice value being 4.52 ]k with proportional lattice constants perpendicular to the surface, and various spacings between the two outermost oxygen layers. This is equivalent to assuming that each oxygen atom, with the exception of the atoms on the surface, is surrounded by four tetrahedral oxygen atoms at identical O O distances. Although, in bulk hexagonal ice, the O O distances for the molecules in the bilayer are slightly different than the O - O distances between bilayers [45], wc have assumed they are identical in the initial model computations. In Section 7.1, we attempt to examine the effect of relaxing this constraint on the best-fit structure. Hydrogen was omitted in most of the LEED calculations, since its contribution is small, partly due to its small scattering cross-section, but primarily in view of the inherent disorder present in its asymmetric positions. We attempted to model

196

N. Materer et aL / Su~Jbce Science 381 (1997) 190 210

the H atom positions by occupying the midpoints of all O O bonds with H atoms having large vibration amplitudes (Debye temperature of 1200K) in order to represent the averaged smeared-out H positions. In addition to hydrogen atoms at the midpoints, calculations were also performed with H in " t o p " sites on the outermost O's of a full bilayer termination to represent surface H atoms not involved in hydrogen bonding. A final calculation was done on a model with H atoms occupying only the top sites to simulate detectable top-site H atoms on the surface with undetectable disordered H atoms in the bulk. These calculations resulted in an Rp-factor value of approximately 0.40. Because these values do not represent an improvement in the Rp-factor over models which do not include hydrogen, we have not attempted any additional hydrogen containing models. This result is not surprising in view of the MD simulations (see Section 6) and suggests disordered hydrogens on the surface as well as in the bulk.

4.2. Analysis of models with a single termination The fit between theoretical and experimental I - V curves for the full- and half-bilayer termination models was quantified by means of the Pendry R-factor (Rp-factor) [47]. The results of structural optimizations as a function of the "bulk" lattice constant d are summarized in Table 1. Three parameters were fit in each case, the inner potential and the perpendicular positions of the first two topmost oxygen atoms in the surface unit cell. A Debye temperature of 300 K (G(Z) of 0.10 A) was used for all O atoms with the exception of the topmost atom in the enhanced vibrational amplitude model. In the case of the full-bilayer model, these two topmost oxygen atoms correspond to Layers 1 and 2 of Fig. 2. For the half-bilayer model, because there is no Layer 1, we fit the perpendicular positions of the oxygen atoms in Layers 2 and 3. Deeper layers were tried but not included in the fit due to the relatively small improvements they made in the resulting Rp-factor. In the model with the enhanced vibrational amplitudes (a(Z) in Layer 1 of 0.28 A), the

Table 1 Rp-factor values for different single termination models"

surface

d (A)

Ih-half

4.33 4.49 4.66 4.99

Ic-full

lc-half

0.10

0.28

0.58 0.47 0.47 0.63

0.56 0.43 0.40 0.45

0.43 0.36 0.35 0.62

lh-full 0.10

0.28

0.52 0.42 0.40 0.40

0.42 0.26 0.27 0.50

0.44 0.28 0.29 0.52

aThe first two interlayer spacings optimized as a function of lattice constant d (parallel to the surface, with proportional lattice changes perpendicular to surface): full- and half-bilayer termination of bulk Ic (cubic) and bulk Ih (hexagonal), with a(Z) (the perpendicular rms vibrational amplitude) of the top oxygen atoms equal to 0.10 or 0.28 A for the full-bilayer model. All other layers have ~r(Z) set to 0.10A. The best Rv-factor values are shown in bold. See Section 7.2 for additional details related to the vibrational amplitudes used in these LEED optimizations. The bulk value of d is 4.52 ~..

topmost atom (Layer 1 ) was held fixed, due to the lack of sensitivity to its position in the Rp-factor0 and the next two oxygen atoms (Layers 2 and 3) were included in the search. These parameters are identical to the parameters fit in the half-bilayer model. Assuming bulk-like vibrations (a(Z) of 0.10 A), our analysis clearly favors the half-bilayer termination and systematically favors the Ih form of ice (Rp-factor near 0.28). The Rp-factor can be improved somewhat with a large vibrational amplitude or o-(Z) of the outermost O layer in the fullbilayer termination of Ih ice (Rp-factor near 0.26). The similarity of Rp-factors for the half-bilayer and the full-bilayer termination with enhanced vibrational amplitude indicates that for these large vibrational amplitudes, the outermost O atoms in the full-bilayer model become essentially undetectable by LEED. However, the presence of these O atoms does have a small influence on the resulting Rp-factor values through forward scattering. The optimization of the perpendicular rms vibrational amplitude or o-(Z) of the top oxygen atoms for the full-bilayer model will be discussed in Section 7.2. The Ih and Ic forms of ice have close Rp-factor values (differences of 0.08 to 0.06), since they differ mainly in the lateral positions of the

N. Materer et a L / Surface Science 381 (1997) 190-210

third and deeper O layers. This structural difference between the two forms is approximately 7 ,~ under the surface, which is comparable to the electron's mean free path. We used Pendry's method [47] to generate the range of Rp-factor values into which all models that are indistinguishable from each other should fall. In doing so, we find that the Rp-factors for the cubic model are just outside this range of values. Another consideration in determining whether we have Ih or Ic ice is the six-fold symmetry of the LEED pattern. The Ic model requires averaging rotated domains in order to produce six-fold symmetrical L E E D patterns, assuming they were formed in the initial growth phase of the film, while similar averaging comes more naturally for the Ih model due to the different orientations of successive terraces. Studies of bulk ice grown on cold metal surfaces have found the Ic form of ice at low temperatures; however, the irreversible transition of the Ic form to the Ih form depends on the thermal history of the sample and is found for a range of temperatures around our growth temperature [45]. For all these reasons, we favor the Ih over the Ic form. 4.3. Analysis of models with two terminations For energetic reasons discussed in Section 5, the half-bilayer termination is likely to be energetically less stable than the full-bilayer termination. One possible way, then, to explain the relatively good agreement between the experimental and theoretical I - V curves for the half-bilayer model is a kinetically limited process during the film growth. This includes condensation on the ice surface or sublimation of water molecules before the sample has had time to cool down to 90 K. For the condensation hypothesis, the diffusion along the surface would have to be too slow, at the low measurement temperature of 90 K, to anneal the remaining half-bilayer into a full-bilayer covering part of the surface, thus, possibly resulting in a mixed termination. The possibility of sublimation before the sample has had time to cool down, seems unlikely from a bond counting perspective. The topmost water molecule in a halfbilayer termination would be bonded only through one hydrogen bond to the next layer; while in a

197

full-bilayer termination, the topmost molecule is bonded to three other molecules. Therefore, it is not clear why the top half of the bilayer would evaporate faster than the exposed bottom half. Instead, one might expect that if sublimation was occurring as the sample cooled, it would help to maintain a full-bilayer termination over most of the surface at all times. T P D studies of a fullbilayer of water molecules on the Ru(0001 ) surface have found that the whole bilayer desorbs at the same temperature [49]. To test a mixed termination model, we have carried out automated T L E E D searches in which the diffracted intensities from both the half- and the full-bilayer termination models were added in several predetermined ratios. This addition of diffracted intensities simulates a surface which contains a mixture of half- and full-bilayer terminated domains. The results of these calculations are illustrated graphically in Fig. 3. In this figure, the left and right endpoints represent surfaces containing only the half- or the full-bilayer termination, respectively. For these endpoints only the three parameters described in Section 4.2 were optimized. Between the two endpoints, a surface which is a mixture of the half- and the full-bilayer termination is simulated. Because the model contains a mixture of two surface terminations, the number of parameters increases to five: the inner 0.45

] 3 parameters

/-

0.40

o.35~'

/

0.30

0.22

/ 5 parameters

i

0.20" ~ - ~ •-,~, I ' 0.00 0.20 0.40

~-0.60

'

'

'

I

0.80

. . . .

1.00

Ratio (Bilayer Area/Total Area) Fig. 3. Results of automated T L E E D searches in which the diffracted intensities from both the half- and the full-bilayer termination models of lh(0001 ) were added in a predetermined ratio (abscissa) in order to simulate a surface which contained a mixture of half- and full-bilayer domains. The curve is a best fit parabola.

198

N. Materer et al. / Sur]ace Science 381 (1997) 190 210

potential and the positions of two oxygen atoms per surface. As one increases the concentration of the full-bilayer termination past 25%, the Rp-factor increases to the value of the full-bilayer termination model. Below 25%, the Rp-factor is approximately equal to the Rp-factor of the halfbilayer termination model. However, one must remember that twice the number of structural parameters were optimized for these mixed models. With these extra degrees of freedom one would expect a significant drop in the Rp-factor as seen in the case of disordered CO on P t ( l l l ) [48]. Thus, the complete lack of a definite minimum between the two terminations strongly supports single termination and not a mixture.

slabs of the ice structure. In Ih ice there are two slightly inequivalent types of H-bonds: those inside each bilayer and those between successive bilayers. We will compare their strengths to each other, to the H-bonds in the water dimer and to those in bulk ice. Since the geometry of the surface bilayers differs from bilayers within a bulk ice crystal, we must perform geometry optimizations to minimize the total energy. We end this section with a brief examination of the potential energy surfaces along selected coordinates in order to obtain an initial approximation of the vibrational properties. 5.1. Lanices with hexagonal symmetry 5.1.1. A single bilayer

4.4. Possible models

Our L E E D computations support the Ih ice phase over the Ic phase. They are consistent with either the half-bilayer termination model with bulk-like vibrational amplitudes or the full-bilayer termination model with enhanced vibrational amplitudes. We have eliminated models which contain mixed termination. The addition of hydrogen does not affect this conclusion. It is thus not possible, with the use of Rp-factors alone, to distinguish between the two termination models. A third possible model, which is also indistinguishable from these models, is a full-bilayer termination model covered by a disordered top layer or film. Previous M D simulations have shown that the topmost layer of ice is disordered (Sv<0.5) at temperatures down to 180 K [5]. We now examine the results of both total energy calculations and M D simulations, in order to distinguish between these models. Then, in Section 7, we can pursue further optimizations of the L E E D results.

5. Total energy calculations

We begin our geometry optimizations using the 631G* basis set to minimize the total energy of a single bilayer. Because the CRYSTAL program utilizes periodicity, we must order the hydrogen atoms in the slab even though these atoms should be arranged in a disordered fashion. In an attempt to estimate the consequences of this localization, we have varied the size of the two-dimensional unit cell to allow for different ordered arrangements. The first arrangement is a (1 x 1) bilayer. In this bilayer all the water molecules have identical orientations in each unit cell (see Fig. 4a) and the two faces of the slab are inequivalent. On one face, H atoms point outward from the slab's surface and on the other face, an oxygen lone pair points outward. The face with external hydrogen

a)

b)

Fig. 4. Top view of two possible arrangements of non-bonded

In this section, we shall examine the energy difference between the half- and the full-bilayer terminations, using total-energy calculations. To distinguish between the two models with different terminations, we need to understand the importance of the hydrogen-bonds (H-bonds) in the

H-atoms in a bilayer formed by the lighter colored water molecules as the top layer and the darker molecules as the bottom. In (a) non H-bonded OH bonds are pointing down into the page, thus, resulting in a (1 x 1) unit arrangement, while in (b) adjacent non H-bonded OH bonds point in opposite directions, forming a (~/3 x 1/3)R30° arrangement. The unit cells are outlined with a dashed line.

N. Materer et al. / Sur[itce Science 381 (1997) 190 210 a t o m s is acidic while the face with external electron pairs is basic. A n o t h e r w a y to l o o k at this situation is t h a t the O H pairs which are n o t H - b o n d e d are o r i e n t e d a l o n g the positive c-axis (the (0001 ) direction) in the (1 x 1) bilayer. T h e o t h e r a r r a n g e m e n t is a ('~/3 x ~ J ) R 3 0 ° bilayer (see Fig. 4b). T h e larger unit cell allows a d j a c e n t n o n H - b o n d e d O H b o n d s to p o i n t in o p p o s i t e directions a l o n g the c-axis, thus reducing d i p o l e d i p o l e repulsion between water molecules. In this model, one third o f the external H a t o m s b e l o n g to one face a n d the r e m a i n i n g two thirds to the o p p o s i t e face. F o r the ( 1 x 1 ) a r r a n g e m e n t , we start o p t i m i z i n g the single bilayer slab by v a r y i n g the O O distance between the w a t e r molecules. By d o i n g so, we vary the lateral lattice c o n s t a n t o f the slab with p r o p o r tional changes between the two oxygen layers. This o p t i m i z a t i o n is very similar, in terms o f geometrical p a r a m e t e r s , to t h a t o f the L E E D lattice c o n s t a n t o p t i m i z a t i o n o f Section 4. A t the m i n i m u m , each oxygen a t o m is at the center o f an i n c o m p l e t e t e t r a h e d r o n f o r m e d b y its three n e a r e s t - n e i g h b o r oxygen a t o m s which are at a distance o f 2.91 ( T a b l e 2: r o w labeled "(1 x 1)"). This distance c o r r e s p o n d s to an O - O interlayer spacing o f 0.97 A. Next, keeping the lateral cell d i m e n s i o n s c o n s t a n t , we allow the v a r i a t i o n o f this interlayer spacing, f o r m i n g a (1 x l ) a r r a n g e d c o m p r e s s e d bilayer. T h e o p t i m a l spacing, 0.592 A instead o f 0.97 A, c o r r e s p o n d s to an O O distance o f 2.81 ( T a b l e 2: row labeled "( 1 x 1 ) c o m p r e s s e d " ) . This r e d u c t i o n o f the O - O distances between the layers is a s s o c i a t e d with an increase o f the cohesive energy within the bilayer. Indeed, the energy per H - b o n d increases from 4.2 to 5.1 kcal tool - 1. This value is c o m p a r a b l e to t h a t for the w a t e r dimer, 5 . 6 k c a l m o 1 - 1 , o r to t h a t o f the bulk, 6. 3 k c a l m o l 1. We n o w t u r n to the (~/3 x X/3)R30 ° a r r a n g e m e n t o f w a t e r molecules in a single bilayer represented in Fig. 4b. In the previous system, the (1 x 1) a r r a n g e m e n t , the H a t o m s on the two faces either all p o i n t o u t w a r d o r all p o i n t inward. Because o f this, the i n t e r a c t i o n between two c o u p l e d (1 x 1) arranged bilayers is asymmetric. F o r the (X/3 x ~/3)R30 ° structure, the larger unit cell allows a d j a c e n t O H g r o u p s to p o i n t in o p p o s i t e directions a n d the two faces o f the bilayer are m o r e similar

199

Table 2 Total energies (in kcal mol 1) gained with respect to the free H20 molecules, calculated with the 631G* basis set for the water dimer, the various bilayer models with hexagonal symmetry and the bulk structure a d (A.)

Hbonds

Energy

Energy/ H-bond

Dimer

2.97*

1

5.6

5.6

Bilayers (lxl) (1 x 1) compressed (V3 x X/3)R30°

2.91" 2.91/2.81"* 2.91

3 3 3

12.7 15.4 15.9

4.2 5.1 5.3

4 Layers (1 x 1) bilayer + bilayer mono + bi+mono

2.91 2.91

7 5

27.0 13.2

3.9 2.6

4 Layers-(~f3 x V3)R30 bilayer + bilayer 2.91 mono + bi+mono 2.91

7 5

39.4 20.0

5.6 4.0

Bulk HOH = 109.5~ HOH=107.1 °

4 4

24.9 25.1

6.2 6.3

2.835* 2.835

aThe values labeled with a * are optimized. The geometry of the oxygen atoms in the bulk is tetrahedral. The value labeled with ** is the optimized O-O distance in a bilayer, keeping the lateral cell dimensions fixed at the value found in the (1 x 1) bilayer optimization. Two different HOH angles were tried for the bulk. When this angle differs from 109.5~ the H atoms do not point toward the oxygen atoms but are off the bond axis. d is the O-O distance in ~.. For the (1 x 1) unit cell "H-bonds" is the total number of H-bonds in a unit cell and for the (~/3 x X/3)R30° unit cell it is the total number of H-bonds in the unit cell divided by 3.

to each o t h e r ( b o t h are a m p h o t e r i c ) when c o m p a r e d to the ( 1 x 1 ) bilayers. This structure is m o r e stable t h a n the (1 x 1) bilayer b u t o n l y by 3.2 k c a l m o l - I In the case o f truly d i s o r d e r e d h y d r o g e n b o n d i n g , the total energy should lie between the (1 x 1 ) a n d (~/3 x ~ J ) R 3 0 c energies; each one represents an extreme c r e a t e d b y the o r i e n t a t i o n o f the O H g r o u p s in the lattice. 5.1.2. Coupling o f the bilayers N o w t h a t we have the g e o m e t r y o f the (1 x l ) a n d the (~/3x~x/3)R30 ° a r r a n g e d bilayers, we e x a m i n e thicker slabs f o r m e d f r o m these two basic units. Presented in Table 2 are the results for f o u r

200

N. Materer et al. / SurJace Science 381 (1997) 190-210

layers using the 631G* basis, arranged either as two full bilayers or as one bilayer flanked on both sides by half bilayers, i.e. monolayers. In both cases, stacks formed using the ( ~ x ~ J ) R 3 0 ° arranged unit cell have a lower energy than the (1 x l) arranged stacks. The most notable result is that the slabs made from coupling two bilayers are much more stable (by 13.8 to 19.4kcalmo1-1) than those where a bilayer is flanked on both sides by half bilayers. The energy resulting from the coupling of the bilayers is weak ( 2 7 . 0 - 2 x 1 2 . 7 or 1.6kcalmo1-1) for the (1 x 1) arrangement, but for the (~/3 x ~3)R30 ~ arrangement it is of the order of the H-bond of the bulk (39.4 2 x 15.9 or 7.6kcalmol-~). This energetic difference illustrates the importance of the orientation of the water molecules in ice.

system. The resulting lattice formed from this cml 1 unit cell is shown in Fig. 5. The cell vectors for the corresponding bulk structure are a = b = 4.88 A, c=4.78 A with angles ~=/~= 118.3 ° and 7=58.6 '. In the following calculations, we have used the 631G** basis set. Results for the bilayer (see Table 3: row labeled "2L") confirm that the H-bond energy of

I

oo

,4

8.51A

5.2. Distorted hexagonal lattices We are interested in the relative energy difference between various slab geometries which differ in their terminations, but not in small and difficult to compute energy differences which give rise to the various ice phases and to the exact bulk geometry. The importance of the H-bonds can be illustrated by keeping all the oxygen atoms at the same positions as in the regular bulk ice (assumed to have perfectly tetrahedral coordination) and placing the H atoms opposite to the neighboring oxygen atoms, thus minimizing H-bonding. From this model one obtains a repulsion of 2.5 kcal mol 1 between the water molecules. The importance of the H-bonds in determining the total energy suggests a geometry in which every O H . . - O bond has the same local geometry as that of the dimer. In the dimer using the 631G** basis set (thus including some polarization to each H atom by the addition of a p-orbital), we calculate the OH bond length to be 0.943 ,~ and the O H . . . O distance to be 2.979 A. The H-bond energy for the dimer is 5.5 kcalmo1-1. To allow the formation of equivalent H-bonds similar to that of the dimer, one must make a small distortion in the layer slab geometry. In particular, the c-axis will no longer be perpendicular to the slab, but must be oriented at an angle of 85.8 °. The resulting system has cml 1 symmetry but remains close to the hexagonal

0.<

Fig. 5. Top view of water molecules positions in a cml 1 model bilayer formed by the lighter colored water molecules as the top layer and the darker molecules as the bottom. The two inequivalent oxygen-oxygen distances (4.88 and 4.78 A) are labeled and the unit cell is outlined by a dashed line. Table 3 Total energies (in kcal m o l - ~) gained with respect to the free H20 molecules, calculated with the 631G** basis set for the water dimer, various models with the distorted hexagonal symmetry derived from the dimer geometry and the cmll bulk model" d (A,)

Hbonds

Energy

Energy/ H-bond

Dimer 2L: bilayer 3L: mono + bilayer 4L: 2 bilayers 6L: 3 bilayers

2.98* 2.98 2.91 2.91 2.91

1 3 5 7 11

5.5 14.1 t6.0 30.1 46.0

5.5 4.7 4.0 4.3 4.2

Bulk cml 1 ice H on bond axis HOH = 109.5

2.900* 2.887*

4 4

25.4 25.2

6.35 6.3

~In slab calculations, these structures have cm I 1 symmetry. For the bulk geometry (labeled as "H on bond axis"), the cell vectors are a = b = 4 . 8 8 ~. and c=4.78 A with angles ~=/~= 118.3 and },=58.6 . The values labeled with a * are optimized. We have not distinguished between the two possible orderings of H-bonds in the slabs due to the similarity of the total energies. d is the O O distance in A,. "H-bonds" is the total number of H-bonds per unit cell.

201

N. Materer et al. / Surj~tce Science 381 (1997) 190 210

4.7 kcal mol =1 within the bilayer is close to that of the dimer ( 5 . 5 k c a l m o 1 - 1 ) or to that of the bulk cml 1 ice (6.35 kcal mo1-1). The energies in Table 3 are results for the dimer, several different bilayer geometries and the bulk c m l l ice. The H atoms are arranged in a (1 x 1) manner. The energies for the four-layer (4L) and the six-layer (6L) systems show that the couplings of the bilayers are weak (30.1 2 x 14.1 or 1 . 9 k c a l m o l 1 and 0.5 x (46.0-3 x14.1) or 1.85 kcal mol - 1, respectively). It is possible, as in the case of the hexagonal bilayers, that this energy could become comparable with that of the H-bond of the dimer or that of the bulk if a (~f3 x X/-J)R30' arrangement of the H atoms is used. The comparison of the 6L energy with the 4 L + 2 L or 3 L + 3 L energies gives the energy for cleavage between bilayers or inside a bilayer. The energy of the cml 1 bulk geometry resulting from the adaptation of the ice geometry to the H-bonds of a dimer is slightly more stable than the regular ice with equivalent O - O distances. However, the energy difference is very small and may disappear when the proton exchanges are considered. Nevertheless, these calculations show the importance of the H-bonds in determining the total energy of the system. 5.3. Vibrations

As the last part of our total energy calculations, we examine the potential energy surfaces along selected coordinates of the water molecules in order to obtain vibrational properties. We compute the total energy as a function of the displacement and fit a second order polynomial to obtain force constants. Assuming a harmonic oscillator model for the vibrations we can then calculate vibrational properties. Both the vibrating water and the remaining part of the ice substrate are assumed to be rotationally and translationally rigid; the H 2 0 reduced mass is taken to be m = 1 8 (vibration against a wall) except for the dimer m = 9 (two punctual masses of 18). The assumed rigidity in this vibrational model is clearly inaccurate, as our M D simulations exhibit (see Section 6). Results are displayed in Table 4 for the 4 layer systems (calculations using the 631G* basis set)

Table 4 Calculated vibrational properties for both the four layer and six layer slabsa

Four layer slab Layer l-(1 x l) Layer2-(lxl) Layer3-(lxl) Layer4 ( l x l ) Layer 1 (X/3x'~/3)R30

kf

Po

PI

rmso rmsl

rmsc

22.3 33.5 34.2 11.9 23.1

0.86 0.91 0.91 0.76 0.86

0.12 0.08 0.08 0.18 0.12

0.09 0.08 0.08 0.10 0.09

0.15 0.14 0.14 0.18 0.15

0.13 0.11 0.11 0.16 0.12

51.7

0.95

0.05 0.07

0.12

0.09

Six layer slab

Layer 1 (1 x l)

"For each layer we give the force constant (kf in kcal (mol ~2) 1) and the population of the ground Po and first vibrational PI states along with the classical rms amplitude (in A,) for three cases: for a half vibrational quantum rmso, for the first vibrational state rms I and for a constant energy of 0.18 kcal mol 1 above the minimum energy rinse.

and for the 6 layer systems (the 631G** basis set is used). From the force constant, we calculate the population (normalized such that the sum over all states is one) of the ground and first vibrational states. We also give the classical rms amplitude for three cases: for a half vibrational quantum (ground state zero point energy), for the first vibrational state and for a constant energy of 0.18 kcal mo1-1 (representing the thermal energy at 90 K ) above the minimum energy. These values are approximately equal to those obtained using a quantum mechanical harmonic oscillator and do not include the anharmonicity of the energy profile. When the vibrating molecule belongs to a surface layer, the potential energy surface is smoother toward the vacuum than towards the bulk, thus increasing the vibrational amplitude when compared to the harmonic potential. The vibrations of molecules in the two inside surface layers are clearly smaller than those on the surface of the slab. In the fourth layer (Layer 4), which has only one H atom oriented toward the slab and involved in an H-bond, the water molecules have the largest vibrational amplitude. The vibrations of the opposite face are smaller; the vibrating molecule involves its two hydrogen atoms in the H-bonds with the slab. We have also considered the vibration of a

202

N. Materer et al. / Sur/hce Science 381 (1997) 190 210

single water molecule at the surface of a (~fJ x ~ J ) R 3 0 ° H-atom arranged four-layer structure. The two other water molecules in the unit cell on the surface remain fixed along with the inner molecules. The water molecule that we have chosen has one H atom pointing outward from the slab into the vacuum. The calculated vibrational amplitude (see Table 4: row labeled "Layer 1-(~f3 x ~/3)R30 °'') is very similar to that obtained for the vibration of the full layer. Table 4 also shows results from similar calculations for the 6 layer structure of cm 11 symmetry. Because this geometry has been constructed to maximize the H-bonds, the amplitude of the molecular vibrations is therefore reduced with respect to the hexagonal slabs. There is no difference between the vibration of species on the two opposite surface layers.

hydrogen bonds extend uninterrupted throughout the liquid [51]. The total energy computations confirm this simple argument for the energetic stability of the full-bilayer termination over the half-bilayer termination. The calculation of the energetic stabilization of full- vs half-bilayer termination gives about 14kcalmo1-1 (total for both sides), which is indeed approximately double the H-bond energy of about 7 to 8 kcal t o o l - 1 In these calculations, we also found that the arrangement of molecular orientations contributes noticeably to the total energy through dipole-dipole interactions, but not to the extent of changing our conclusions. Calculations of the potential profile for individually displaced molecular layers also imply enhanced surface vibrations. The possibility of enhanced surface vibrations is considered in more detail in the next section using MD simulations.

5.4. Conclusions from total energy computations 6. Molecular dynamics simulations In this section, we have examined the two possible surface terminations discussed in Section 4, using total energy calculations. One would expect the ice surface to maximize the number of hydrogen bonds. Such maximization of course occurs within bulk ice [45], and also appears to occur at water surfaces and water/molecule interfaces [50]. Qualitatively, one would expect an energy difference between the half- and full-bilayer termination to be equal to the energy of one hydrogen bond. When the bulk lattice is cleaved such that both surfaces terminate in the fullbilayer, one hydrogen bond per surface unit cell is broken, as can be seen by examining Fig. 2. However, if the lattice is cleaved such that both surfaces terminate in the half-bilayer, three hydrogen bonds are broken. Thus, this simple computation gives an energy difference of one hydrogen bond per free surface unit cell per surface in favor of the full-bilayer termination. Our calculations also show that hydrogen bonding is the main factor in the determination of the surface structure. This is true even within liquid water where, on average, each molecule makes 2 to 3 hydrogen bonds with neighboring molecules, forming a random network, such that strings of

Given the Hartree-Fock results, we now examine the results of the molecular dynamics (MD) simulations described in Section 3.3. With these simulations we can obtain dynamic information about the water molecules at the ice surface. For both the full- and half-bilayer termination models, we examine the rms vibrational amplitudes of the surface molecules perpendicular and parallel to the surface. Using the translational order parameter Sv, we address the long-range order of the topmost water molecules.

6.1. Full-bilayer termination Starting with the full-bilayer termination of the lh(0001) surface, we collect in Table 5 the rms amplitudes a(Z) of motion perpendicular to the surface, translational order parameters (Sv) and rotational order parameters (SR) for the upper eight layers calculated using the MD method. The most important result of this table is that a(Z) of the first layer is much larger than that of the other layers and also much larger than the bulk value of the zero-point uniaxial rms vibration amplitude (0.13 A at 90 K ) [45]. Thus, the MD simulations

N. Materer et al. / SutJiwe Science 381 (1997) 190-210 Table 5 M D simulation results of the surface with a full-bilayer termination for the rms amplitudes of motion of the O-atoms perpendicular a(Z) and parallel or(X) (equal to a(Y)) to the surface a Layer

a(Z )

G(X )

ST

SR

1 2 3 4 5 6 7 8

0.24 0.20 0.19 0.16 0.16 0.15 0.15 0.14

0.26 0.26 0.24 0.24 0.21 0.21 0.18 0.18

0.90 0.9 l 0.92 0.92 0.94 0.94 0.95 0.95

0.66 0.86 0.87 0.88 0.89 0.90 0.90 0.91

a6 is given in A. for the first 8 layers, or the first four bilayers. The translational order parameter ST and the rotational order parameter SR are also shown.

lend support to the assumption that, in the full bilayer termination, the water molecules in the upper layer of the first bilayer are essentially "invisible" in the LEED experiments (in the sense of Table 1 and Fig. 13), due to their enhanced amplitude of motion perpendicular to the surface. In addition, the amplitudes of motion parallel to the surface in the first layer is much enhanced over that of the lower layers. The a(Z) value calculated for the first layer is actually not quite large enough to make the upper layer molecules totally "invisible" in the LEED experiment: some residual sensitivity to the outer molecules should remain, as implied in Section 4. Indeed, close examination of Table 1 and Fig. 13 shows that the upper layer is not completely "invisible". In addition, the MD result for a(Z) is likely to be a lower bound for a(Z) of the "experimental ice", for several reasons. In the M D simulations, the contribution of the intramolecular vibrations to the oxygen amplitudes is absent (the water molecules are treated as rigid rotors). Furthermore, the lattice constant parallel to the surface which fits the experimental results best is enhanced over the bulk value, indicating the presence of some strain in the surface: possibly the ice layer is not yet thick enough for the surface molecules to be unaware of the presence of the metal surface underneath. In the presence of such strain, one may expect the hydrogen bonds to be

203

somewhat weakened, resulting in increased amplitudes of motion. According to the MD simulations, the rms amplitude a(Z) slowly decreases when going into the bulk (see Table 5). Note that the amplitude in layer 2 is similar to that in layer 3, and the same is true for layers 4 and 5, and so on. The amplitudes of motion perpendicular to the surface are more or less the same for the molecules of the lower layer or half-bilayer within a given bilayer as those for the molecules in the upper layer or half-bilayer of the bilayer below. Because the hydrogen bonds between these molecules (for example the molecules in layers 2 and 3 of Fig. 2) are oriented perpendicular to the surface, the motions of the molecules perpendicular to the surface are highly correlated in these half-bilayers. This correlation would reduce the effective rms amplitudes as seen by LEED electron scattering computation, because correlated motions between two atoms are not taken into account. This could possibly explain why LEED sees a much faster decay of the rms vibrational amplitudes (see Section 7.2) than the MD simulations. Finally, the amplitude of motion as calculated classically may be somewhat too small due to the neglect of zero-point vibrational motion. MD simulations that incorporate quantum effects in the vibrations of the surface molecules [18] show that quantum effects may make a small, but significant contribution to the amplitude of motion of the molecules in the outermost layers. The enhanced vibrational amplitude in the upper layers is also visible in Fig. 6, where we show a density profile along the Z-direction (see also Section 3.3). As is clearly evident, the two layers in the first bilayer have started to merge. From a dynamical perspective, it will no longer be possible to make the assignment of a given molecule to one of the layers within the first bilayer on the basis of the instantaneous Z-coordinate of the molecule only. In contrast, the layers within the third, fourth and fifth bilayer are basically distinct. The molecules in the sixth bilayer are clearly affected by the presence of the rigid layers underneath, as is evident by the sharp peak produced by the lower layer. Turning our attention to the translational and

204

N. Materer et al. /Surjace Science 381 (1997) 190 210 1.0

10 z

D >0£

< o/

.

.

.

i

.

.

.

.

i

.

.

.

.

0.8

8

I

.

.

.

.

LAYER 1

>0:1

surface

~<

6

0.6

I-

I--

•(

4

Eo

2

v >I--

.

v >-

Z LJ E]

Z Ld rm

0.4

0.2 0.0

o 0

15 20 z (ANGSTROM)

25

Fig. 6. The density profile of a surface with a full-bilayer termination (surface is at right) is shown along the direction perpendicular to the surface at T=90 K. The atomic density is in arbitrary units, and the distance is in A,. Only the moving layers are shown.

•

1.0

30

.

-0.5

0.0

0.5

.0

1.0 i-

LAYER 2

0.8 >or"

0.6

rotational order parameters (see also Section 3.3 presented in Table 5, we find that the surface o f ice is m u c h more ordered at 90 K than it is at 1 9 0 K and higher [5]. The translational order parameters for the first and second layers are 0.90 and 0.91 at 90 K, respectively; at 190 K they are in the range o f 0.4 to 0.6 (see Ref.[5]). The translational order parameter m a y be seen as a structure factor, and the large value found for the first and second layers indicates that these layers should appear ordered in a L E E D experiment. The M D results for the surface with a full-bilayer termination are consistent with the full-bilayer termination f o u n d by the L E E D I - V curves analysis, the first layer being "invisible" due to the large vibrational amplitudes o f the surface molecules and not due to disorder. At 90 K, the rotational order parameter SR is 0.66 for the first layer; it is in the range o f 0.2 to 0.3 at 190 K. Nevertheless, SR is m u c h lower for the first layer (0.66) than it is for the second and deeper layers (0.86 and higher), indicating that the extent o f rotational disorder is m u c h larger in the first layer. The rotational disordering is visualized in Figs. 7 and 8, where we plot distributions ("densities") o f cos 0 and ¢ (see Section 3.3. As can be seen from Fig. 7, the distribution o f cos 0 a b o u t the ideal lattice values (0.5 and - 0 . 5 ) is

I---

<

0.4

0.2

~

0.0

/ .0

-0.5

0.0 COS THETA

0.5

1.0

Fig. 7. The resulting distribution of cos 0 (in arbitrary units) for the water molecules in the first two layers at T=90 K. 0 is the polar angle the dipole vector of water makes with the surface normal. These results are obtained for a surface with a fullbilayer termination. much wider for the first layer than it is for the second. Likewise, the C-distribution, shown in Fig. 8, is also broader for the first layer. The increased amplitude o f m o t i o n a ( Z ) and the enhanced rotational disorder in the first layer should be a consequence o f the disruption o f the hydrogen bonding network at the surface [17]. In bulk ice, the molecules are firmly held in place by a strongly directional tetrad o f hydrogen bonds, which is disrupted at the surface. In the ideal (0001) surface with full bilayer termination, the molecules in the outermost layer only participate in three hydrogen bonds. As a consequence, they m a y try to move inwards to form an extra

205

N. Materer et al. / Surface Science 381 (1997) 190 210

layer for the present experiments enhanced value o f ~r(Z).

1.0 i-

LAYER

0.8

I

J

6.2. Half-bilayer termination

>.r~

a:: < 0.6

1

I-

O.4 >0.2 Z Ld C2~

0.0 1

1

z

. 0

........

0.8

i . . . . . . . . .

2

i . . . . . . . . .

,3

, . . . . . . . . .

4

, . . . . . . . . .

5

, . . . . . . . . .

6

i

LAYER 2

>r,F"

0.6

2

0.4 yF--

N

is still the

0.2

Z LLI

0.0 0

1

\~. 2 3 4 PHI (RADIANS)

5

6

Fig. 8. The resulting distribution of ~b (in arbitrary units) for the water molecules in the first two layers at T=90 K. ~ is the azimuthal angle (in radians) which defines the projection of the dipole vector of water on the surface. These results are obtained for a surface with a full-bilayer termination.

h y d r o g e n bond. This m a y lead to strains and some disruption o f the hydrogen bonding network. At T = 190 K, this strain also leads to high rotational mobility o f the molecules, as is evident from the small value calculated for the N M R correlation time (~NMR of the order o f 100 ps [5]). We f o u n d it impossible to calculate a converged value for rNMR (it should be larger than 1 ns) at a temperature o f 90 K. Nevertheless, the low value o f SR for the first layer and the wide angular distributions indicate an enhanced rotational mobility o f the molecules in the first layer. O f course, the most important consequence o f the strain in the surface

To investigate whether the resulting experimental L E E D I - V curves could also possibly be consistent with the half-bilayer termination model, we performed additional M D simulations o f such a surface. The vibrational amplitudes of m o t i o n and order parameters obtained from the calculations are summarized in Table 6. We find similar amplitudes as in the full-bilayer termination (see Table 5). However, the most important result of this simulation is that the translational order parameter for the first layer is n o w very small, approximately equal to one divided by the n u m b e r o f molecules in the first layer. This amounts to a zero structure factor for the positions o f the oxygen atoms relative to one another, meaning that this layer should not be seen as a crystalline layer in the L E E D experiment, but as a m o r p h o u s or diffusing. In the simulation o f the surface with a full-bilayer termination, the distance between the lowest (first visible) layer in the uppermost bilayer and the highest layer in the bilayer below (the distance between layers 2 and 3) is 2 . 7 2 A , in reasonable agreement with the L E E D results. In the surface with the half-bilayer, the distance between layer 1 and the uppermost layer in the bilayer below has collapsed to a value of 2.37 A. Table 6 MD simulation results of a surface with a half-bilayer termination for the rms amplitudes of motion of the O-atoms perpendicular a(Z) and parallel or(X) (with is equal to a(Y)) to the surfacea Layer

a(Z)

a(X)

ST

SR

1 2 3 4 5 6 7

0.22 0.20 0.16 0.16 0.15 0.15 0.14

0.26 0.18 0.18 0.17 0.17 0.16 0.16

0.03 0.86 0.88 0.92 0.92 0.95 0.95

0.00 0.51 0.85 0.86 0.86 0.90 0.90

act is given in A for the first 7 layers, i.e. one half-bilayer plus the next three bilayers. The translational order parameter ST and the rotational order parameter SR are also shown.

206

N. Materer et al. / SmJbce Science 381 (1997) 190 210

This short distance is well outside of the range of values (2.80_+0.10 A) determined by the L E E D L V curve analysis. Thus, we conclude that our M D results for the surface with the half-bilayer t e r m i n a t i o n are inconsistent with the interpretation of the L E E D ~ V curves as being due to a surface with a half-bilayer t e r m i n a t i o n . A plot of the density profile of the surface with half-bilayer t e r m i n a t i o n is given in Fig. 9. As can be seen, a low-density overlayer is present (at 29 A,), which was formed by one molecule j u m p ing out of the u p p e r layer. As can also be seen, the bilayer below the first layer is much more disordered t h a n the second bilayer in the simulation of the surface with a full-bilayer t e r m i n a t i o n (see Fig. 6), the two layers within this bilayer being almost no longer distinct in the half-bilayer case. As can be seen from Table 6, the rotational order parameters for the u p p e r m o s t layer in the surface with a bilayer t e r m i n a t i o n are zero, indicating that this layer, in addition to being t r a n s l a t i o n ally disordered, is also rotationally disordered a n d is n o t at all similar to an ice like layer. A more qualitative view is provided in Figs. 10 a n d 11, where we plot distributions ("densities") of cos 0 a n d ¢ (see S e c t i o n 3 . 3 ) . As carl be seen from Fig. 10, the d i s t r i b u t i o n of cos 0 is almost flat for the u p p e r m o s t layer, as it should be for a liquid-

z

t--

z

_

LAYER 1

0.8

< 0.6 02 F--

0.4

0.2 Z Ld C)

0.0

.0

-0.5

0.0 COS THETA

0.5

.0

Fig. 10. The resulting distribution of cos 0 (in arbitrary units) for the water molecules in the first two layers at T=90 K. 0 is the polar angle the dipole vector of water makes with the surface normal. These results are obtained for a surface with a halfbilayer termination. 1.0t. ........ , ......... , ......... , ......... , ......... , . . . . . . . , z

0.8 I

LAYER 1

>02

0.6 F-

<~ 0.4 )-

0.2 0.0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 1 2 .5 4 5 6 PHI (RADIANS)

8

>Od

Fig. 11. The resulting distribution of ~ (in arbitrary units) for the water molecules in the first two layers at T-90 K. ~ is the azimuthal angle (in radians) which defines the projection of the dipole vector of water on the surface. These results are obtained for a surface with a half-bilayer termination.

surface

< 02

6

<

4

IF

2 Z LIJ d:)

. . . .

:>{M

Z LIJ K3

10

q

1.0 (Y? F-

0 0

5 20 Z (ANGSTROM)

25

50

Fig. 9. The density profile of a surface with a half-bilayer termination (surface is at right) is shown along the direction perpendicular to the surface at T=90 K. The atomic density is in arbitrary units, and the distance is in A. Only the moving layers are shown.

like layer. O n the other hand, the d i s t r i b u t i o n of ~b is certainly not flat, indicating that some orientations are favored t h r o u g h the interaction with the ice-like bilayer below.

7. R e f i n e m e n t s t o t h e L E E D

structure

After c o m b i n i n g the L E E D results with the total energy a n d molecular dynamics results, we con-

N. Materer et al. / Surface Science 381 (1997) 190-210

clude that the Ih(0001) model, with full-bilayer termination and enhanced vibrational amplitudes for the top oxygen atoms is clearly the best model for our ice film. The total energy calculations and the MD simulations strongly suggest that the apparent halfbilayer termination (the invisibility of the upper half-bilayer in LEED) is due to large amplitude vibrational motion of the uppermost molecules and not to disordering of the top layer. The half-bilayer termination is ruled out by total energy arguments and by the MD results, which show that the half bilayer should be invisible in LEED due to disordering. In this section, we refine this initial model of Section 4 to obtain the bestfit structural parameters.

7.1. Refinement of the lattice parameters In bulk ice, the O O distances within the bilayer are slightly different from those between bilayers [45]. In the LEED computations (see Table 1 in Section 4) we have assumed that all O - O distances are identical. We now attempt to examine the effect of relaxing this artificial constraint on the best-fit structure. For the best-fit Ih(0001) model, with full-bilayer termination and enhanced vibrational amplitudes, we have explored different lateral lattice constants from 4.33 to 4.99 A parallel to the surface (see Table 1), the bulk ice value being 4.52 A, with proportional lattice constants perpendicular to the surface, and various spacings between the two outermost oxygen layers. This procedure resulted in an optimal lateral lattice constant of 4.62 ,~. By allowing the lateral lattice constant to change proportionally to the lattice constant perpendicular to the surface, we have enhanced the sensitivity of the Rp-factor to changes in its value. Because fitting does not take into account that the O-O distances for the molecules in the bilayer are slightly different from those between bilayers in bulk ice, we have also performed LEED calculations for models in which we varied only the lateral lattice constant with the perpendicular lattice constant equal to the above value. Fig. 12 shows the resulting Rp-factor variation as the lateral lattice constant took on values between 4.33 and 4.99 A with a minimum at 4.62 (the bulk value being 4.52 A). This optimization

207

0.55 ~1~

\

0"45-1

0.40~

~ 0.35~ ¢~ 0.30 ~ 0.25 ~ 0.20 ~ 4.20

? "k

/

--~

. T

"x,,,,,x ~ . . ! 4.40

~

/ .

I 4.60

l 4.80

5.00

2D Lattice Constant Fig. 12. Rp-factor variation of lh(0001 ) with the 2D lattice constant. The lateral lattice constant was varied from 4.33 to 4.99 parallel to the surface with the lattice constant perpendicular to the surface fixed at its optimal value (the bulk value being 4.52 A). For each lattice constant value (marked by an "x"), the outermost two interlayer spacings were optimized. The solid line is a best-fit parabola.

scheme has the effect of varying both the in-plane O-O bond distance and bond angle. Unfortunately, the Rp-factor has low sensitivity to further refinements of the lattice constants, due to the relatively small energy range of the LEED data. For the final structure, we have assumed that all bulk O O bond distances are identical.

7.2. Enhanced vibrations In Section 4.2, we used a bulk rms vibrational amplitude of 0.10A and an enhanced amplitude of 0.28 A without justification. Now that the fullbilayer termination model with enhanced vibrational amplitudes is supported by both the total energy computation and the MD simulations, we can examine these parameters in more detail. Fig. 13 shows the Rp-factor variation vs the inverse rms perpendicular amplitude of the topmost oxygen. Plotting the Rp-factor variation vs the inverse rms permits a better graphical display. The value of zero for the inverse rms amplitude corresponds to the half-bilayer termination model. The best-fit value (rms amplitude of 0.25 to 0.33 A) is close within the range of values predicted by the MD simulations of Section 6. However, the error bar on this determination is relatively large. The Rp-factor was relatively insensitive to the vibra-

208

N. Materer et al. / Sur[ace Science 381 (1997) 190 210

0.45

Table 7 aOxygen coordinates (in A) obtained by LEED for one 2D unit

0.40-~ 0.35

cell of the best-fit full-bilayer termination with enhanced vibrational amplitudes for the top oxygen atoms model of ice lh(0001 )

o.3oi

Atoms

X

O1 02 03 04

0.00 1.33 1.33 0.00

0.00 2.3l 2.31 0.00

0.00 0.99 3.79 4.80

Repeating bulk atoms 05 06 07 08

0.00 1.33 1.33 0.00

0.00 2.3l 2.31 0.00

7.63 8.57 11.40 12.34

0.00

0.00

7.55

....

4.00 4.00

2.31 2.31

0.00 0.00

--

0.25 0.20~ 0.00

~ 2.00

T 4.00

r 6.00

8.00

10.00 12.00

Y

Z (insensitive) _+0.05 _+0.06

l / rms z

Fig. 13. The Rp-factor as a function of the inverse rms perpendicular amplitude of the topmost water layer. The bulk value of l/rms is about 10. The value of zero for the inverse rms amplitude corresponds to the half-bilayer termination model and the minimum value is close within the range of values predicted by the MD simulations. The data points are marked with an "x" and the curve is a best fit parabola.

--

Bulk repeat vector

2D lattice vectors

tional amplitudes o f deeper layers. It is possible that the correlated m o t i o n o f the 2nd layer with the 3rd layer (see Fig. 2) described in Section 6 reduces the effective rms amplitudes seen by the L E E D electron. The L E E D c o m p u t a t i o n s do not take correlated motion into account. We have chosen 0 . 1 0 A for the bulk rms perpendicular amplitude because it corresponds to experimental values for bulk Ih ice [45].

aX and Y are coordinates parallel to the surface (repeating with the given 2D lattice vectors) while Z is perpendicular to the surface. The top 4 bilayers are explicitly listed, atoms Ol through 0 4 belonging to layers 1 through 4 of Fig. 2. The 3rd and 4th bilayers (05 0 8 ) , together with the bulk repeat vector, define the bulk structure. One unit cell is given by its 2D lattice vectors along with its repeat vector for the bulk. Due to the enhanced vibrational amplitudes, the Rp-factor was insensitive to variations in the position of the top oxygen atom; thus, its position was not optimized.

7.3. The f i n a l structure

For the best-fit full-bilayer termination with enhanced vibrational amplitudes we obtain an Rp-factor o f 0.26. The lattice constant parallel to the surface, with proportional lattice changes perpendicular to the surface, was varied to obtain an optimal value o f 4 . 6 2 ± 0 . 1 0 A , This corresponds to an ideal hexagonal unit cell with edges a = 4 . 6 2 and c = 7 . 5 6 A,, within the ice film, c o m p a r e d with bulk values o f 4.52 and 7.37 A, respectively. Thus, the lattice o f the ice film m a y be slightly expanded by 2_+3%, possibly due to the mismatch between the ice film parallel to the surface and the Pt substrate; a match would require a 6.4% lateral expansion of the ice. The coordinates o f the final best-fit structure are given in Table 7; the corresponding best fit I - V curves are shown in Fig. 14.

8. Conclusions We find the ice surface at 90 K to be terminated by strongly vibrating molecules in the top half o f the outermost bilayer, such that the molecules would be undetectable by L E E D . This model results in the best agreement between the theoretical and experimental intensity vs energy ( I - V ) curves and is consistent with total energy considerations based on maximizing the n u m b e r o f hydrogen bonds. The model is also consistent with M D simulations, which show large amplitude vibrational m o t i o n o f the uppermost water molecules in an otherwise stable and well ordered bilayer termination, and bulk-like interlayer spacings.

A( Materer et al. /' Surlctce Science 381 (1997) 190-210

[,0] & [O,l]

..~'~ ~

[1,11 ...

-

.~

[2,01 & 10,21

Experiment Theory lT ~ I 30 50

~ T

~ '

70

7

~

•

90

r~-,--,- 1- : ~ 7 • 1 110 130 150

Energy Fig. 14. The best-fit intensity vs energy ( ~ V ) curves for the Ih(0001) full-bilayer termination model with the strongly vibrating (or otherwise disordered) molecules in the top half of the bilayer. The coordinates of the final best-fit structure, which yields an Rp-factor of 0.26, are given in Table 7.

Acknowledgements U.S. gratefully acknowledges financial support by the Deutsche Forschungsgemeinschaft (DFG). G.J.K. acknowledges support by the Koninklijke Nederlandse Academie van Wetenschappen (KNAW). This work was supported in part by the Director, Office of Energy Research, Office of Basic Energy Sciences, Materials Sciences Division of the US Department of Energy under Contract No. DE-AC03-76SF00098. References [1] H.H.G. Jellinek, J. Colloid Interface Sci. 25 (1967) 192. [2] J.G. Dash, Contemp. Phys. 30 (1989) 89. [3] F.P. Bowben, D. Tabor, in: The Friction and Lubrication of Solids, Clarendon Press, Oxford, 1964, part II, ch. VII. [4] M. Elbaum, M. Schick, Phys. Rev. Lett. 66 (1991) 1713. [5] G.-J. Kroes, Surf. Sci. 275 (1992) 365. [6] A. Goto, K. Akiya, T. Hondoh, Y. Furukawa, T. Shimura, I. Takahashi, J. Harada, J. Crystal Growth 129 (1993) 491. [7] M. Elbaum, S.G. Lipson, J.G. Dash, J. Crystal Growth 129 (1993) 491.

209

[8] A. Lied, H. Dosch, J.H. Bilgram, Phys. Rev. Lett. 72 (1994) 3554. [9] H. Dosch, A. Lied, J.H. Bilgram, Surf. Sci. 327 (1995) 145. [10] J.G. Davy, G.A. Somorjai, J. Phys. Chem. 55 (1971) 3624. [11] V. Buch, J.P. Devlin, J. Phys. Chem. 98 (1993) 4195. [12] J.D. Graham, J.T. Roberts, J. Phys. Chem. 98 (1994) 5974. [13] B. Rowland, M. Fisher, J.P. Devlin, J. Phys. Chem. 97 (1993) 2485. [14] P.A. Thiel, T.E. Madey, Surf. Sci. Rep. 7 (1987) 211. [15] H. Ogasawara, J. Yoshinobu, M. Kawai, Chem. Phys. Lett. 231 (1994) 188. [16] G. Held, D. Menzel, Surf. Sci. 316 (1994) 92. [17] T.A. Weber, F.H. Stillinger, J. Phys. Chem. 87 (1983) 1983. [18] E. Batista, H. J6nsson, to be published. [19] X. Xia, M.L. Berkowitz, Phys. Rev. Lett. 74 (1995)3193. [20] X. Xia, L. Perera, U. Essman, M.L. Berkowitz, Surf. Sci. 335 (1995) 401. [21] U. Starke, N. Materer, A. Barbieri, R. D611, K. Heinz, M.A. Van Hove, G.A. Somorjai, Surf. Sci. 287 (1993) 432. [22] U. Starke, K. Heinz, N. Materer, A. Wander, M. Michl, R. D011, M.A. Van Hove, G.A. Somorjai, J. Vac. Sci. Technol. A 10 (1992) 2521. [23] D.F. Ogletree, G.S. Blackman, R.Q. Hwang, U. Starke, G.A. Somorjai, J.E. Katz, Rev. Sci. lnstrum. 63 (1991) 104. [24] M.A. Van Hove, W. Moritz, H. Over, P.J. Rous, A. Wander, A. Barbieri, N. Materer, U. Starke, D. Jentz, J.M. Powers, G. Held, G.A. Somorjai, Surf. Sci. Rep. 19 (1993) 191 : "Barbieri-Van Hove Automated Tensor LEED Package" available from M.A. Van Hove. [25] R. Dovesi, C. Pisani, C. Roetti, M. Cause, V.R. Saunders, CRYSTAL 88, QCPE program No. 577, Bloomington, IN, 1989. [26] P.C. Hariharan, J.A. Pople, Theoret. Chim. Acta. (Bed.) 28 (1973) 213. [27] W.J. Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 56 (1972) 2257. [28] W.J. Hehre, R.F. Stewart, J.A. Pople, J. Chem. Phys. 51 (1969) 2567. [29] J.S. Binkley, R.A. Whiteside, R. Kishnan, R. Seeger, H.B. Schlegel, D.J. Defrees, J.A. Pople, Gaussian 80, QCPE No. 406, Bloomington, IN. [30] S.J. Szalewicz, W. Cole, W. Kolos, R.J. Bartlett, J. Chem. Phys. 89 (1988) 3662. [31] J.C. White, E.R. Davidson, J. Chem. Phys. 93 (1990) 8029. [32] L. Ojam~te, K. Hermansson, R. Dovesi, C. Roetti, V. Saunders, J. Chem. Phys. 100 (1994) 2128. [33] J.M. Besson, P. Pruzan, S. Klotz, G. Hamel, B. Silvi, R.J. Nelmes, J.S. Loveday, R.M. Wilson, S. Hull, Phys. Rev. B 49 (1994) 12540. [34] E. Cota, W.G. Hoover, J. Chem. Phys. 67 (1977) 3839. [35] J.D. Bernal, R.H. Fowler, J. Chem. Phys. 1 (1933) 515. [36] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. lmpey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [37] W.L. Jorgensen, J.D. Madura, Mol. Phys. 56 (1985) 1381. [38] O.A. Karim, A.D.J. Haymet, J. Chem. Phys. 89 (1988) 6889.

210 [39] [40] [41] [42] [43]

N. Materer et al. / Surface Science 381 (1997) 190-210

S.W. de Leeuw, J.W. Perram, Mol. Phys. 37 (1979) 1313. E.R. Smith, Mol. Phys. 57 (1986) 793. D.E. Parry, Surf. Sci. 49 (1975)433. D.E. Parry, Surf. Sci. 54 (1976) 195. H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [44] D. Fincham, Mol. Simul. 8 (1992) 165. [45] D.S. Eisenberg, W. Kauzmann, The Structure and Properties of Water, Oxford University Press, New York, 1969.

[46] E. Walley, in: P. Schuster, G. Zundel, C. Sandorfy (Eds.) The Hydrogen Bond, North-Holland, Amsterdam, 1976, p. 1425. [47] J.B. Pendry, J. Phys. C 13 (1980) 937. [48] G.S. Blackman, M.-L. Xu, D.F. Ogletree, M.A. Van Hove, G.A. Somorjai, Phys. Rev. Lett. 61 (1988)2353. [49] G. Pirug, C. Ritke, H.P. Bonzel, Surf. Sci. 241 ( 199l ) 289. [50] Y.R. Shen, private communication. [51] F.H. Stillinger, Science 209 (1980) 451.