J . Phys. Chem. 1989, 93, 6171-6180
6171
Computer Simulation of the Water/Platinum Interface E. Spohrt Max-Planck-Institut fur Chemie (Otto-Hahn-Institut). P.O.B. 3060, 0 - 6 5 0 0 Mainz, FRG (Received: January 18, 1989; In Final Form: March 17, 1989)
A molecular model for the water/platinum interface has been devised. It includes the surface corrugation of the metal and orientationally anisotropic water-metal interactions obtained from quantum chemical cluster calculations. Barriers for the surface diffusion and for the reorientation of a single water molecule on the quadratic (100) face of the face-centered-cubic platinum crystal are discussed. The flexible BoppJancs6-Heinzinger water model describes the water-water interactions, and the platinum-platinum interactions are described by a single force constant. Molecular dynamics simulations of a water lamina confined by (100) platinum surfaces have been performed using these interaction potentials. The structure is discussed on the basis of one-particle density profiles and solvent pair correlation functions. The surface-induced structural inhomogeneity ranges up to distances of IO A. In the center of the lamina water properties are bulklike. Hydrogen bonding in the vicinity of the interface is only slightly reduced relative to the bulk. The orientational structure of water is strongly influenced by water-water interactions and is considerably different from the preferential orientation according to the water-platinum interaction potential. It leads to a dipolar potential drop across the interface of 1 . 1 V.
1. Introduction An understanding of the properties of liquid water in contact with metal surfaces on a molecular basis is of eminent interest in both electrochemistry and catalysis. Recent experiments on water adsorption on clean single-crystal surfaces of transition metals under ultrahigh vacuum have provided valuable data on the structure, especially the orientational structure, and on the dynamical properties like hindered translational, rotational, and vibrational motions of adsorbed water. The concentrations of water range from nearly isolated molecules a t submonolayer coverages up to icelike solid.'-8 Molecular dynamics (MD) simulations are able to reproduce the general behavior of the work function change upon water adsorption at around monolayer coverages and thus allow the interpretation of this behavior on the basis of experimentally inaccessible distribution functions? This work is extended now to the liquid water/platinum two-phase system which represents more the standard electrochemical situation than the single crystal with adsorbed water under vacuum. M D data allow the estimation of the range and magnitude of changes induced by the presence of the metal on structural and dynamical properties of water using a realistic description of both water-water and water-surface interactions. Analytical theory as a means of the description of the water/metal interface is presently elaborated at the level of the dipolar hard sphere or the Stockmayer model (Lennard-Jones point dipole) for the The metal is usually described by the jellium model. While this approach is adequate for s-p metals it fails to describe any of the characteristic binding properties of water on the catalytically and electrochemically more interesting d metals like e.g. platinum.'* Systematic computer simulation studies have also been performed for the Stockmayer liquid in contact with hard wall^.'^-^^ Recent simulation work has been concerned with a regular array of adsorption sites on an otherwise hard wall,Is and the effect of an orientationally anisotropic potential (which is realistic for the water-metal interactions) on the density distribution close to a surface has been treated by integral equation theory and Monte Carlo simulations.I6 The point dipole model is, however, not adequate for a realistic description of water" since it does not account for the specific anisotropic hydrogen-bond properties which are included in the various point charge models for liquid water. Studies with the Stockmayer model are thus not capable of a realistic description of the competition and balance between adsorption forces due to the metal and the water-water interactions. Many computer simulation studies of liquid/gas or liquid/solid systems have been performed using realistic water models together
+
'Present address: Department of Chemistry, University of California, Irvine, CA 92717.
0022-3654/89/2093-6171$01.50/0
with a simplified description of the water-surface interactions. In most cases the solid was modeled as a molecularly smooth wall. A review of this work can be found in ref 18. In the present approach a more realistic model for water-platinum interactions is introduced and discussed, and results obtained by means of Molecular dynamics (MD) simulations for a water lamina in contact with the (100) lattice plane of the face-centered-cubic platinum crystal are presented. The flexible Bopp-Jancsb Heinzinger (BJH) water modelI9 has been used since it was shown in extensive simulations on ionic hydration that it adequately describes the effect of intermolecular forces on water properties (for recent reviews see ref 20 and 21). The flexibility of the water molecule, although computationally much more expensive than a rigid model, was regarded to be a necessary constituent of a realistic model in order to investigate possible changes in intramolecular geometry upon adsorption as suggested by quantum chemical calculations.22 Furthermore, the comparison with vibrational modes from electron energy loss spectroscopic measurements is rendered possible. Platinum was chosen for several reasons: (i) A variety of experimental data is a ~ a i l a b l e . ~ ~ ~ ~ ~ ~ - ~
(1) Crowell, J. E.; Chen, J. G.; Hercules, D. M.; Yates, Jr., J. T. J. Chem. Phys. 1987,86, 5804. (2) Bange, K.; Grider, D.; Sass, J. K. Surf. Sei. 1983, 126, 437. (3) Peebles, D.E.; White, J. M. Surf. Sei. 1984, 144, 512. (4) Langenbach, E.; Spitzer, A.; Liith, H. Surf.Sei. 1984, 147, 179. (5) Fusy, J.; Ducros, R. Surf. Sei. 1986, 176, 157. (6) Nobl, C.;Benndorf, C. Surf. Sei. 1987, 182, 499. (7) Bange, K.; Madey, T. E.; Sass, J. K.; Stuve, E. M. SurJ Sei. 1983, 183, 334. (8) Thiel, P. A.; Madey, T. E. Surf. Sei. Rep. 1987, 7, 211. (9) Spohr, E.;Heinzinger, K. Ber. Bunsen-Ges. Phys. Chem. 1988, 92, 1358. (10) Badiali, J. P. Mol. Phys. 1985, 55, 939. (1 1) Badiali, J. P.; Boudh-Hir, M. E.; Russier, V. J . Electroanal. Chem. 1985, 184, 41. (12) Kornyshev, A. A. In The Chemical Physics of Solvation. Part C. Studies in Physical and Theoretical Chemistry; Dogonadze, R. R., Kilmin, E., Kornyshev, A. A., Eds.; Elsevier: Amsterdam, 1988; Vol. 38, p 61. (13) Lee, S. H.; Rasaiah, J. C.; Hubbard, J. B. J. Chem. Phys. 1986,85, 5232. (14) Lee, S.H.; Rasaiah, J. C.; Hubbard, J. B. J . Chem. Phys. 1987,86, 2383. (15) Caillol, J. M.; Levesque, D.; Weis, J. J. J . Chem. Phys. 1987, 87, 6150. (16) Russier, V.; Rosinberg, M. L.; Badiali, J. P.; Levesque, D.; Weis, J. J. J . Chem. Phys. 1987,87, 5012. (17) Goldman, S.;Backx, P. J . Chem. Phys. 1986,84, 1986. (18) Spohr, E.;Heinzinger, K. Electrochim. Acta 1988, 33, 1211. (19) Bopp, P.; Jancd, G.; Heinzinger, K. Chem. Phys. Lett. 1983,98, 129. (20) Heinzinger, K. Physic0 1985, 131B, 196. (21) Bopp, P. Pure Appl. Chem. 1987, 59, 1071. (22) Bauschlicher, Jr., C. W. J . Chem. Phys. 1985, 83, 3129.
0 1989 American Chemical Society
6172 The Journal of Physical Chemistry, Vol. 93, No. 16, 198'9
Spohr convenience an ansatz of pairwise additive interactions between oxygen and platinum and hydrogen and platinum has been made. Based on extended Hiickel molecular orbital calculation^^^ of a water molecule on top of a five-atom platinum cluster (several geometries have been investigated), the pair interaction potential has been adjusted to a set of exponential functions. The analytical form of the water-platinum pair potential is given by dH20-Pt = d'O-Pt(rO-Pt,PO-Pt) + dH-Pt(rH,-Pt) + d'H-Pt(rH2-Pt) (1)
2
i
d
6 AzIA
Figure 1. Total water-metal interaction potential as a function of the oxygen-surface distance Az for water on top of a platinum atom (full line), on a bridge site (dashed), and on a &fold hollow site (dotted). The dipole moment vector of the water molecule points away from the surface in all three cases.
(ii) Dissociation of water, which can presently not be treated adequately in a computer simulation study, can be-at least at low temperatures-neglected.26 (iii) Surface corrugation and preferential bonding through the oxygen atom were modeled on the basis of barriers for diffusion on the surface and reorientation of the water molecules taken from extended Huckel calculation^?^ The interaction potentials employed in the simulation are discussed in section 2. Section 3 gives the technical details of the simulation. Questions concerning the range of the effect of water-metal interactions on water properties are addressed in section 4. The effect on the liquid structure is discussed in section 5 on the basis of pair correlation functions which are compared to those of bulk water, and the changes in hydrogen bonding are investigated in section 6. The orientational order of the water molecules relative to the metal surface is discussed in section 7 and leads to an idealized structural model of the interfacial region based on several prevailing water configurations. This is followed by a short summary of structural results in section 8. A future communication will be concerned with dynamical properties represented by autocorrelation functions and their Fourier transforms. Comparison with spectroscopicdata of water adsorbed in various concentrations on clean single-crystal platinum surfaces will be made in order to assess reliability and realism of the employed microscopical model.
2. Interaction Potentials Water-water interactions are described by a modified central force potential (BJH p ~ t e n t i a l )which '~ has been tested in simulations of bulk water, clusters, and electrolyte solutions. This model describes water structure and dynamics and especially the effect of intermolecular interactions on the intramolecular geometry and vibrations fairly well and thus lends itself to a study of the effect of an extended planar metal surface on water properties. The platinum-platinum interactions are modeled by a harmonic nearest-neighbor potential2*in order to account for the coupling between motions of the water molecules and lattice vibrations. The development of a water-metal interaction potential function was guided by the requirements that it should (for an isolated water molecule on the extended surface) firstly correctly describe the adsorption site on top of a metal atom which is found in all quantum chemical calculations of water on metal surfaces known to the author. Secondly, it should exhibit a minimum for orientations with the oxygen atom closer to the surface than the hydrogen atoms. The latter feature is consistent with both the negative work function change observed on transition-metal surfaces upon water adsorption2-' and the theoretical studies of the bonding of water to metal ~ l u s t e r s . ~ UFor ~ ,computational ~~~ (23) Sexton, B. A. Surf. Sci. 1980, 94, 435. (24) Fisher, G. B.; Gland, J. L. Surf. Sci. 1980, 94, 446. (25) Wagner, F. T.; Moylan, T.E. Surf. Sci. 1987, 191, 121. (26) Ibach, H.; Lehwald, S . Surf. Sci. 1980, 91, 187. (27) Holloway, S.; Bennemann, K. H. Surf. Sci. 1980, 101, 327. (28) Black, J. E.; Bopp, P. Surf. Sci. 1984, 140, 275. (29) Bauschlicher, Jr., C. W. Chem. Phys. Lefr. 1987, 142, 71.
where = [1894.2 exp(-I.l004r) 1886.3 exp(-l,0966r)lf(p) + lo6 exp(-5,3568r)[ 1 - f ( p ) ] (2) and dH-Pt= 1.7 142 exp(-l.2777r) (3) with
= exp(-0.5208~~) (4) The energy is given in units of joules, and r and p are given in angstroms. The sitesite distances are denoted by r. p is the length of the projection of the distance vector onto the surface plane. The pair potential daPt is attractive for larger Pt-0 distances only if p is small, leading to the preferential adsorption site on top of a platinum atom with an adsorption energy of 35.7 kJ/mol on the (100) surface (see below), in agreement with the extended Hiickel calculation^.^^ Calculations of water interactions with Al,31Ni," and Cu30 clusters using different methods and cluster sizes yield binding energies of the same order of magnitude, namely, 51.1, 33.7, and 36.6 kJ/mol, respectively. These adsorption energies agree reasonably well with estimates obtained from thermal desorption spectra which range from 40 to -70 kJ/mol for most transition metal^.^"^.^^ Due to the p dependence of the potential, those platinum atoms with small overlap between oxygen and platinum orbitals ( p large) repel1 the water molecule. This is consistent with the finding by Ribarski et al.30 that the binding energy decreases approximately by a factor of 2 in going from the H20-Cu to the H20-Cu5 cluster. A quantitative fit to the data in ref 27 was not attempted as the Hiickel method has severe limitations and can be taken only as a guideline. Specifically, it is expected to be most reliable in the region of the potential minimum as for very small distances exchange interactions are not properly accounted for and for large distances the finite cluster size comes into effect. In the following the total energy Vof one water molecule above the (100) surface of an extended Pt crystal will be discussed. Figure 1 shows the variation of the total water-metal interaction energy as a function of the distance Az from the crystal surface. (Throughout this paper the surface is defined as the equilibrium position of the top layer of platinum atoms.) Curves have been calculated separately for water on top of a surface atom (full line), on a bridge site between two platinum atoms (dashed), and on the Cfold hollow sites of the (100) face (dotted) with an orientation where the dipole moment of the water molecule points away from the surface which is the energetically most favorable orientation (see below). In agreement with the quantum chemical calculations, the minimum of the total energy is found for the on-top site. For bridge and hollow sites the position of the energy minimum shifts considerably to larger distances while its depth (-35.7, -16.3, and -1 2.1 kJ/mol for on-top, bridge, and hollow sites, respectively) agrees qualitatively with the values of -51.1, -24.1, and -14.5 kJ/mol calculated by Muller and Harris" for water on top of AI(100). Holloway and B e n n e n ~ a n ndo ~ ~not give energies for f(p)
(30) Ribarski, M. W.; Luedtke, W. D.; Landman, U. Physica 1985,832, 1430. (31) Muller, J . E.; Harris, J. Phys. Rev. Left. 1984, 53, 2493. (32) Anderson, A. B. Surf. Sci. 1981, 105, 159. (33) Estiu, G.;Maluendes, S . A,; Castro, E. A,; Arvia, A. J. J. Phys. Chem. 1988, 92, 2512. (34) Madey, T. E.; Yates, Jr., J. T. Chem. Phys. Left. 1977, 51, 77. (35) Benndorf, C.; Nobl, C.; Thieme, F. Surf. Sci. 1982, 121, 249.
Computer Simulation of the Water/Platinum Interface
i
i
6
8 Azla
Figure 2. Total water-metal interaction potential as a function of the oxygen-surface distance Az for different orientations of water on top of
a platinum atom: dipole vector pointing away from the surface (full line) and pointing toward the surface (dotted); dipole vector in the surface plane and proton-proton vector parallel to the surface (dashed) and perpendicular to the surface (dash-dotted). the bridge site. The quantum chemical calculations indicate that the transition state for surface diffusion is on a bridge site which is, therefore, modeled correctly by the interaction potential model employed here. However, no experimental evidence is known to the author which proves this prediction. In Figure 2 the Az dependence of the total water-metal interaction energy is depicted for four orientations of water on top of a platinum atom. The energy minimum is found for an orientation with the dipole moment vector pointing away from the surface (full line). The most unfavorable orientation is the one with the dipole moment pointing toward the surface (dotted). The exponentially repulsive hydrogen-platinum interactions lead to a negligible contribution to the surface corrugation. This is in accordance with the quantum chemical result that the rotation of water about the A 1 4 axis is practically unhindered on an A19 cluster.” Figure 3a-d shows the surface corrugation in such a way that the minimum total potential energy (obtained from the variation with the distance Az for a given pair (x,y))as a function of the coordinates x and y parallel to the surface for fixed orientations equivalent to the ones in Figure 2 is plotted. It can be seen that the surface roughness is the more pronounced the deeper the minimum of the potential is (see Figure 2). This is again in accordance with the results of Muller and Harris on A1(100), namely, that the variations of the orientational energy are relatively weaker on the hollow and bridge sites than on the on-top sites. Thus, the energy difference between molecules on different sites with the same orientation becomes smaller. The barriers for surface diffusion along the (1 10) direction ( x = y ) are considerably larger than the thermal energy kT and amount to 19.4, 14.8, 12.3, and 9.1 kJ/mol for a, b, c, and d, respectively.
3. Details of the Simulations
w
The rectan ular basic box has side lengths L, = L,, = 19.6 8, and L, = 45 and is shown in Figure 4b. The interval -12.7 C z C 12.7 8, is occupied by 305 BJH water molecules, and in the regions 12.7 IIzI I22.5 A 550 platinum atoms are located. As a consequence of the introduction of periodic boundary conditions in all directions, the platinum subsystem can be regarded as one infinitely extended slab with a thickness of 19.6 8,which interacts with water molecules on both surfaces. The values of Lx and L,, have been chosen in accordance with the periodicity of the face-centered-cubic platinum lattice with a lattice constant a = 3.92 8,. The crystal slab consists of 5 X 5 X 5 face-centered-cubic unit cells, each containing four atoms and the (100) plane being parallel to the (xy) plane of the periodic box. In Figure 4a the surface geometry is sketched by circles and crosses denoting surface atoms and second layer atoms, respectively. The system was prepared in such a way that the water density in the inner region of the lamina (-5 C z C 5 8,)is 1 g/cm3. After an extensive equilibration period of 15 ps the simulation extended over 15000 time steps of 0.50 fs with occasional temperature scaling and subsequently of 35 000 steps of 0.25 fs without rescaling the velocities. This is equivalent to a total elapsed time of 16.25 ps. Structural properties have been obtained as averages
The Journal of Physical Chemistry, Vol. 93, No. 16, 1989 6173 over the whole simulation period whereas dynamical properties were calculated only from the final 35 000 steps. No differences in structural properties between the two intervals have been found. The average temperature was 349 K. The high temperature has k n chosen in order to obtain a sufficient statistical accuracy and not to trap the system in a localized region of phase space. The Ewald summation method for two-dimensional system^'^-^^ has been used to calculate the long-range Coulombic interactions. Details are given in the Appendix. For the non-Coulombic interactions the shifted force method39 was used with a cutoff distance of Lx/2.
4. Density Profiles and Dipolar Potential Drop In previous simulations of water laminas it was usually found on the basis of density fluctuations that a distance of the order of 20 8, is sufficient to decouple the interfaces (see e.g. ref 18 and references therein). Figure 5 shows the results for the relative oxygen and hydrogen atom density profiles p(’)(Az). A pronounced and well-separated first peak in p g ) at Az = 2.5 8, reflects the directly adsorbed water overlayer, A shoulder on the long distance side of the peak indicates that a fraction of the water molecules is pushed out of their energetically most favorable adsorption sites (see also section 7). This feature has been found even more pronounced for a water m ~ n o l a y e r .The ~ reason for this behavior is that although the Pt-Pt nearest-neighbor distance (2.77 A) is very similar to an oxygen-oxygen distance of 2.76 and 2.85 8, in ice and liquid water, respectively, the quadratic arrangement of the adsorption sites is incompatible with a water molecule forming two linear hydrogen bonds at a time. Therefore, an energetically more favorable situation obviously occurs when some of the water molecules are slightly pushed out of the overlayer. A second peak in the oxygen density profile occurs at Az = 5.4 A. It is less pronounced and seems to be produced by those water molecules which are hydrogen bonded to the first layer. For even larger distances the density is practically constant except for statistical noice and equal to the density of bulk water. For nonpolar liquids a more pronounced layered structure has been found than for water (see e.g. ref 40). The first maximum of the hydrogen atom density profiles coincides with that of the oxygen profile but is less pronounced than the latter. This indicates a strong preference for an orientation of water molecules parallel to the surfaces. A low but well-pronounced maximum at Az = 4.4 A coincides with the minimum of the oxygen profile. It is followed by a broad maximum at Az = 5.9 A. The maximum at 4.4 A is due to hydrogen bonds between first and second layer (see below) whereas the maximum at 5.9 8, mainly follows the increase in the oxygen density at 5.5
A.
The high barriers for surface diffusion of water (see Figures 1 and 3) lead to a localization of water molecules on the preferred on-top sites. The range of this lateral ordering can best be investigated by an expansion of the oxygen atom density po(x,y,Az) into a Fourier series according to
po(x,y,Az) =
,,kpij(Az)
,,=-m
COS
(y ) (7) COS
(5)
1/(L,LY6z)S L X ’ h xSLy’2dySi\r+az’2dz~ p o ( ~ , y , z 3x -1J2
-LJ2
AZ-62’2
cos
(y) (y ) cos
(6)
(36) Heyes, D. M.; Barber, M.; Clarke, J. H. R. J . Chem. SOC., Faraday Trans. 2 1977, 73, 1485. (37) Heyes, D. M.; van Swol, F. J . Chem. Phys. 1981, 75, 5051. (38) Heyes, D. M. J . Chem. Phys. 1983, 79, 4010. (39) Streett, W. B.; Tildesley, D. J.; Saville, G. In Computer Modeling of Matter; Lykos, P.,Ed.; American Chemical Society: Washington,DC, 1978; ACS Symp. Ser. No. 86, p 144. (40) Lee, C. Y.; McCammon, J. A.; Rossky, P. J. J. Chem. Phys. 1984, 80, 4448.
6174 The Journal of Physical Chemistry, Vol. 93, No. 16, 1989
Spohr
Figure 3. Total water-metal interaction potential V(x,y;Az=Az&) for different orientations of water on top of a platinum atom: dipole vector pointing away from the surface (a) and pointing toward the surface (d); dipole vector in the surface plane and proton-proton vector parallel to the surface (b) and perpendicular to the surface (c). Azmin= Azmin(x,y)is the energetically most favorable oxygen surface distance for a given position (x,y).
a)
t'
1011I
f
--0.4
-0.51
1
,
1
x/v
1.5 -
V/kJmol"
- 30
'1 l o I
0.5
-30
\ +'
L
J
1
1
1
2
4
6
0
1
IO MA
1
1
1
'
1
1
2
4
6
0
II IO d z l a '
Figure 5. Oxygen (a) and hydrogen (b) density profiles p ( ' ) , normalized
Fourier coefficient pss/poo (c) (for definition see text), dipole moment density pg (d), running potential drop x (e), and average water-metal interaction energy ( f ) (together with the energetically most favorable curve (dashed) as shown in Figures 1 and 2). Az is the distance from the surface.
r~~
Figure 4. (a) Sketch of the arrangement of surface platinum atoms
(circles) and second-layer platinum atoms (crosses) in the xy plane of the basic cell. The dashed lines are the periodic boundaries of the basic cell at x,y = f9.8 A. The [OlO] and [OI 11 directions of the crystal are indicated. (b) Sketch of the basic tetragonal simulation cell. The water molecules are located in the center of the box, and the platinum atoms are represented by the circles. The coordinate system which is used throughout the paper is sketched in the center. where L, and Ly are the dimensions of the elementary cells, 62 = 0.1 A, and i, j are integers. In this denotation poo is the oxygen density profile in Figure 5a. As the box parameters L, and Ly have been chosen to be 5 times the lattice constant a of the crystal,
contributions to p(x,y,Az) from those pij where i, j are multiples of 5 are a measure of lateral order. Figure 5c shows the coefficient ps5/poowhich equals -1 if the oxygen atom is located on an on-top site while it is +1 on the 4-fold hollow sites. The division by poo eliminates the effect of density variation perpendicular to the surface. Throughout the regime of the first peak of p g ) (Figure sa) the normalized Fourier coefficient is nearly constant. The large negative value of =-0.8 suggests a strong correlation between oxygen atoms and platinum equilibrium sites. Between the first and second layer the Fourier coefficient increases to zero and then drops again to a value of ~ - 0 . 2for the second layer. This shows that the water molecules in the second layer are more evenly distributed over the ( x y ) plane than those of the first layer. However, the negative sign of the coefficient indicates that some water molecules are located just on top of those in the first layer. This is in keeping with the small intermediate maximum in pfi' at 4.4 8, and the distance of 2.9 A between the first and second layer of oxygen atoms (Figure 5a). Beyond the second layer the
Computer Simulation of the Water/Platinum Interface
The Journal of Physical Chemistry, Vol, 93, No. 16, 1989 6175
coefficient is zero within the limits of statistical accuracy. This shows that there is no remaining lateral correlation with the surface. Higher order Fourier coefficients support this interpretation and are not shown here. The variation of the dipole moment component perpendicular to the surface (Figure 5d) determines the electrical properties of the interface. The solvent contribution to the interfacial potential drop can be calculated from the dipole density p,, through
with co the dielectric permittivity of the vacuum. It has been demonstrated that with the present model good agreement can be achieved with experimental work function changes for partially filled water mono- and bilayer^.^ It is evident from Figure 5d that only the first layer of water molecules exhibits a significant net dipole moment. The dipole density is, however, smaller by about a factor of 10 than it would be for a completely ordered layer at this density. This is in keeping with the coincidence of the first maxima of the oxygen and hydrogen density profiles which suggests an orientation of the molecular plane nearly parallel to the surface (see above). It demonstrates that water-water interactions strongly determine the behavior of the adsorbate layer as the most favorable orientation of the water molecule would be perpendicular. In and beyond the second layer the dipole density is small and scatters significantly. There is, however, considerable anisotropy in the dipole moment distribution, as will be discussed below. Figure 5e shows the running integral ~(Az). In accordance with the dipole density, it can be seen that only the first layer contributes to the potential drop whose asymptotic value is 1.1 f 0.2 v. The surface potential x ( m ) is larger than the value of 0.8 V found for a water coverage 8 = 1.4 (Le., 1.4 water molecules per surface platinum atom).9 The reason for this is mostly because in the case of 8 = 1.4 the second layer contributes negatively to x due to a preferentially antiparallel dipole alignment relative to the first layer which is due to hydrogen-bonding interactions. The effect of the additional bulk water in the present simulation is to suppress this behavior. There is a large cancellation of contributions due to different distances. Because of the relatively short simulation time, the dipole density does, however, fluctuate significantly. The high surface potential is consistent with the general behavior of the work function change upon water adsorption”’ on clean metal surfaces where the multilayer adsorption (comparable to the system studied here) usually leads to an additional decrease of the work function by about 20% from the level a t approximately monolayer coverage. Figure 5f supports the conclusions drawn above concerning the importance of water-water interactions on the basis of the simulation-averaged water-metal interaction energy. The observed interaction energies are typically about 8 kJ/mol (-3kT) higher than that of the energetically most favorable configuration (dashed; see also Figure 1 ) . This is in agreement with the observation that water in the first layer is significantly tilted from the most stable configuration. It is seen that even with the strongly anisotropic water-metal interactions employed in the present work most static one-particle properties decay within three layers (“8 A) from the wall, which corresponds roughly to the range of the water-metal interaction potential. The simulation box seems thus to be large enough to simulate two isolated interfaces.
5. Pair Correlation Functions The strong water-surface interaction together with the corrugation of the surface gives rise to much more drastic changes in water structure than the ones observed in previous computer simulations of water near smooth unpolar surfaces (see references in ref 18). The arrangement of the adsorption sites in a quadratic grid with a spacing of 2.77 A allows on the one hand the formation of hydrogen bonds between pairs of adsorbed molecules but is on the other hand not compatible with the tendency to form a local tetrahedral network like in bulk water (see also section 6). Hence,
IJ+- -
1. \
I
,
l
I
I
,
I
I
with p(‘) the local i-particle density and 7, and F2 the Cartesian coordinates of the two particles. In Figure 6 the radial distribution functions (RDFs) goo(r), gOH(r),and gHH(r) are shown for reference particles i in the adsorbed layer (full line) and in the second layer (dashed). The radial structure around the molecules in the second layer is very similar to that in the third layer and in the bulk liquid at 336 K!l The layer of adsorbed water molecules (full curve) shows a behavior significantly different from the bulk but similar to the one found in a previous simulation42with different water-water and water-metal interaction potentials. The first-neighbor peak is strongly enhanced, but its position is not shifted from its bulk value of 2.85 A. The crystal matrix obviously supports the oxygenoxygen arrangement in spite of the fact that the platinum-platinum distance in the crystal is slightly shorter than that. Beyond the first-neighbor shell several new maxima occur in goo at about 3.9, 5.6, 6.3, and 8.7 A together with a shoulder around 8 A while the maximum found in bulk water around 4.5 A completely vanishes. The new maxima are induced by the periodicity of the surface adsorption sites which are positioned at 2112, 2, 5’12, 8II2, 3, and times the nearest-neighbor distance. Due to the subdivision of the RDFs into layers, the relation g, = gee does not generally hold. Here, however, the oxygen-hydrogen RDFs goH and gHO (the former are the RDFs around an oxygen in a specified subsystem, the latter the RDFs around a hydrogen) are identical. Therefore, only goH is shown in Figure 6. Only minor differences between all three subsystems are found for the first-neighbor peak at about 2 A. The second maximum is shifted by more than 0.1 %, to smaller distances and seems to consist of two overlapping peaks. This effect is due to the quadratic arrangement of the adsorbed molecules that leads with high probability to hydrogen bonds directed along the two surface diagonals x = y and x = y ([Oll] and [Oli] direction; cf. Figure 4a). While in a tetrahedral arrangement the distance between oxygen and second hydrogen is about 3 . 3 A, the corresponding (41) Jancsd, G.; Bopp, P.; Heinzinger, K. Chem. Phys. 1984, 85, 377. (42) Spohr, E.; Heinzinger, K. Chem. Phys. Lett. 1986, 123, 218.
6176
The Journal of Physical Chemistry, Vol. 93, No. 16, 1989
Spohr
Figure 7. Anisotropic oxygen-oxygen pair correlation functions goo(p,z) for (a) the adsorbate molecules (2.1 < Az, I4.2 A), (b) the molecules in the second layer from the surface (4.2 C Azl I6.4 A), and (c) the molecules in the center of the lamina (6.4 A < Azl). p is the transversal and z the normal part of the interatomic distance, and Az, is the distance of the reference atom from the surface.
distance on the quadratic lattice is about 3.0 8,. The observed value of the position of the second peak in the adsorbate layer lies between these bounds. In the intermediate layer the R D F is again similar to the bulk. The first intermolecular peak of gHH in the adsorbate layer includes, for the same reason as in goH, a large number of short distances. Contrary to the second layer and the bulk, there exists also a shoulder at -2.8-3.0 8, which is due to linear O-H--O-H arrangements along the [Ol11 or the [Oli]direction on the surface. The far-reaching small-amplitude oscillations of goH and gHH can be better understood together with Figure 8 (see below). Figure 7 demonstrates the degree of anisotropy of the pair correlation functions for water molecules close to the wall (a), in the intermediate layer (b), and in the center of the lamina (c). It shows goo(p,z) which is a representation of the P C F as a function of two components of the distance vector. p = roo sin t9 is the component of the interparticle distance parallel to the surface, and z = roo cos 8 is the component perpendicular to the surface. 8 is the angle between the surface normal vector ( z direction) and the 0-0 vector roo. Owing to the fact that in a 2-dimensional representation of a P C F fewer events are counted per histogram bin than in a 1-dimensional one, and that the statistical uncertainty is thus increased, the data had to be smoothed before plotting. As a consequence, the heights, especially of the first-neighbor peaks, are significantly reduced relative to Figure 6 . The peak positions, however, remain unaffected. From Figure 7a one concludes that the change of g m ( r ) in the adsorbate layer beyond the first-neighbor shell is due solely to correlations parallel to the surface, which show the characteristic far-reaching oscillations. In addition, there is an increased probability of finding a neighbor in the direction away from the surface as indicated by the peak at z = 3 8, and p 0 A. This observation is reflected by the increased probability for second-layer molecules to find a neighbor at z = -3 8,and p 0 8, (Figure 7b). These oxygen atoms are responsible for the negative sign of ps5 in the second layer. Furthermore, the intermediate maximum of the hydrogen atom density profile at 4.4 8, is obviously due to hydrogen bonds between these molecules and the first layer. In the center of the lamina the oxygen-oxygen PCF is isotropic (Figure 7c). The observed dependency of the correlation from the angle between the 0-0 vector and surface normal is negligible within the limits of the statistical uncertainty (as especially for p 0 only very few atom pairs contribute to g(p,z)). In order to get a closer view on the correlations between adsorbed molecules, PCFs ga8(x,y)are plotted in Figure 8. x and y are the components of the interparticle vector along the two laboratory frame axes parallel to the surface (the [OOl]and [OlO] directions of the crystal, respectively; cf. Figure 4), and both atoms are located in the adsorbed layer. The far-ranging symmetrical arrangement of high maxima in the oxygen-oxygen P C F indicates the formation of a quadratic overlayer commensurate with the crystal periodicity which was already inferred from Figures 6 and 7. The 0-0 correlation is obviously not due to water-water interactions but a consequence of the quadratic arrangement of adsorption sites. Due to the water-surface potential employed (in which the hydrogen atoms do not experience major barriers on the surface), the hydrogen atoms experience the surface corrugation only through the oxygen-crystal interactions. Consequently, goH shows
-
-
Figure 8. Adsorbate-adsorbate pair correlation functions gaB(x,y)with a,0 = 0, H and 2.1 < Azl, Az2 5 4.2. x and y are the projections of the interatomic distance on the x and y directions of the laboratory coordinate system. Az, and Az2 are the distances of the atoms from the
surface. a much smaller degree of correlation. There are groups of four maxima in goH, each corresponding to one of the four different possible directions of an intramolecular 0-H bond along the two diagonals ( x = y , x = -y) relative to an oxygen atom at the origin, a consequence of the arrangement of adsorption sites. The hydrogen-hydrogen PCF exhibits only very little correlations. The hydrogen atoms are disordered in a similar manner as in proton disordered ice. 6. Hydrogen Bonding It was pointed out in section 4 that the lateral interactions between the adsorbed molecules lead to a decrease in adsorption energy relative to the values for isolated molecules. On the other hand, the strong changes in the PCFs indicate the perturbation of the hydrogen-bond network in the interfacial region. It is therefore appropriate to study the hydrogen-bonding properties in more detail, as they are, by means of their effect on the molecular orientation, essential for the understanding of electrochemical properties. In the following two water molecules will be termed "hydrogen bonded" if the 0-0 distance is smaller than 3.35 A and if the angle t9,, between the intramolecular 0 - H vector and the intermolecular 0-0 vector is smaller than 20'. This geometrical hydrogen-bond definition has been applied
Computer Simulation of the Water/Platinum Interface I
I
I”
6
2
The Journal of Physical Chemistry, Vol. 93, No. 16. 1989 6177
4
6
8
10 Azl%
Figure 9. Number of neighbors nN (r I3.35 A) and number of hydrogen-bonded neighbors nH (for definition see text) as a function of the distance Az from the surface.
-10
20
60
60
80
100
a,
Figure 10. Distribution of hydrogen-bond angles OHB (the angle between the intramolecular oxygen-hydrogen vector and the oxygen-oxygen vector for 0-0 distances of O. This is in agreement with the overall slightly negative dipole density in this range (Figure 5d). Beyond 6.4 A there is another inversion of p ( m 190H).All other anisotropies that persist beyond Az = 6.4 A are only weakly pronounced. The broad distributions do not allow a clear distinction of a prevailing configuration in this range and are therefore not discussed in detail. They show, however, the decay of orientational order with increasing distance from the metal surface. Lee et aL40 presented a detailed analysis of the orientational structure of an ST2 water lamina between smooth (9-3) Lennard-Jones walls. They showed that the orientation close to the wall is consistent with a thermally disturbed ice I structure where the c axis is perpendicular to the surface. In the present case the picture is more complicated because of the strong orienting forces exerted by the platinum lattice. A pronounced change of all distributions with distance is observed within the first water layer. The changes in the distribution of cos 9, are not easy to understand. Prevailing orientations can, however, be uniquely characterized by monitoring the probability of a water molecule to form simultaneously angles between the surface normal and the molecular dipole vector 9, and proton-proton vector 9, with the surface normal. These probabilities are depicted in the form of 3-dimensional plots in Figure 16 for the first three distance slabs in Figures 12-1 5. Close to the surface, where the distributions have only one maximum, a single configuration prevails with 9, ;= 75' and 9, ;= 90' (top). Subsequently, with increasing distance from the surface 9, becomes larger (center) and is on the long distance side of the adsorbed water layer (bottom), =120°. These two configurations agree with an icelike structure as in ref 40. Simultaneously, a second kind of molecule with 9, = 60' and 9, from 0 to 30' from the normal develops and contributes positively
to the surface potential drop x (Figure 5e). Figure 17 sketches the four major idealized configurations found in Figure 16. The most probable configuration I generates the maximum at 9, ;= 90' and 9, = 75'. Configurations I1 and 111 contribute to the maximum at 19, ;= 0-30' and 9, ;= 60' (center and bottom), and IV corresponds to the maximum at 9, ;= 90' and 6, ;= 120'. Obviously, hydrogen bonding between any two configurations is possible. However, on the quadratic grid of the platinum surface a continuous hydrogen-bond network cannot be formed easily. This is consistent with the decrease of the number of hydrogen bonds close to the surface (see Figure 9). It is interesting to note that configuration I is one in which the water molecule is bonded to the surface roughly in the lone-pair direction in agreement with the minimum-energy orientations obtained for isolated water molecules on Cu and AI c l ~ s t e r s , ~ ~ ~ ~ ~ although neither the water-water (contrary e.g. to the S T 2 modela) nor the water-wall interactions explicitly make use of a lone-pair site. Obviously, the water-water interactions favor this orientation as the water-platinum interactions alone would lead to a preference of the antidipole orientation. There is some controversy in the literature about the range of the orientational structure in computer simulations of water laminas. The results presented here clearly indicate that at a distance beyond about 10 A from a corrugated surface with anisotropic watersurface interactions all orientational distributions are isotropic. Lee et aL40 also find isotropic distributions at distances of about 10 A using a much less strongly interacting surface model. Sonnenschein and H e i n ~ i n g e as r ~ well ~ as Valleau and Gardner& found that a lamina consisting of 6 or even 12 layers is not thick enough to effectively decouple the two interfaces. They found a prevailing probability for orientations with the dipole moment parallel to the surface at all distances. As isotropic distribution functions in the center of the lamina have been found for simulations with the interactions potentials used in this work in which the shifted force method has been used to truncate the long-range interaction^,“^ it is not likely that the various different treatments of long-range forces are responsible for the differences. It seems more likely that the length of the simulation and the choice of initial configurations are more important.
8. Summary and Conclusions The water-crystal interaction potential employed in the present study describes the behavior of a single water molecule on a platinum surface reasonably well. The adsorption energy is within the range of experimentally determined adsorption energies of water on transition-metal surfaces (between -40 and -70 kJ/ mo15J4934235).The on-top adsorption agrees with quantum chemical calculation^,^^-^^^^^^^ and the binding occurs through the oxygen atoms in agreement with both quantum chemical calculations and work function measurement^.^-^ The metal-induced changes in water structure range up to a distance of about 10 8,which is the same as the range of the direct water-crystal interactions. On the contrary, simple liquids show in general pronounced density oscillations which range, due to packing effects, significantly farther than the direct interactions. Pair correlation functions and orientational distribution functions (44) (45) (46) (47)
Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1974, 60, 1545. Sonnenschein, R.; Heinzinger, K. Chem. Phys. Lett. 1983, 102, 550. Valleau, J. P.; Gardner, A. A. J . Chem. Phys. 1987, 86, 4162. Spohr, E. Unpublished results.
6180
The Journal of Physical Chemistry, Vol. 93, No. 16, 1989
show that the center of the lamina is isotropic, and the radial distribution functions agree with those obtained from a bulk M D simulation with the same water model. The reason that the changes in water structure are only short ranged may be attributed to the relative strength of the hydrogen bonds between water molecules compared with water-platinum interactions. It was shown above that water is able to form a quadratic adsorbate overlayer on the Pt( 100) surface with only a small disturbance of the number and stability of hydrogen bonds. The water-water interactions also determine the orientation for molecules in the surface layer and thus the surface potential. It is interesting to note that water molecules close to the metal surface are bonded roughly in the lone-pair direction, in spite of the fact that the water-crystal interaction potential favors an antidipole orientation. Consequently, the knowledge of the exact binding geometry of a single water molecule on metal surfaces (by means of different quantum chemical calculations both "electrostatic" antidipole and "covalent" lone-pair bonding configurations have been f o ~ n d ~ ' *is~ not * ) crucial for the description of an extended water phase in contact with a metal surface. The orientation of water molecules is not uniform across the adsorbate layer. With increasing distance from the surface different preferential orientations occur, each, however, with a broad distribution around the average. Hydrogen bonding between these molecules is possible in agreement with the small reduction of the number of hydrogen bonds close to the surface. Furthermore, there is considerable hydrogen bonding between the adsorbate layer and subsequent layers. The present work demonstrates that an understanding of electrochemical properties of waterltransition metal interfaces on a microscopical basis requires the knowledge of both realistic water-water interactions (to describe hydrogen bonding) and water-metal interactions which can be taken from quantum chemical cluster calculations. The orientational ordering can only be understood as a complex interplay between these two kinds of interactions which can be conveniently revealed by molecular dynamics simulations.
Acknowledgment. It is a pleasure to thank K. Heinzinger and P. Bopp for many helpful discussions. Financial support by Deutsche Forschungsgemeinschaft is gratefully acknowledged.
Appendix In order to reduce the C P U time for the calculation of the Coulombic interactions, the formulas for the Ewald method in the 2-dimensional case (see Appendix of ref 36) have been evaluated for q = 0.2L, and then fited to the following heuristic function: p3 v =fi PI - + f2 + -
p2
p4
Spohr
fi and& are functions depending on the reduced distance r = r'/L; 1
fl(r) = - + 0.041r - 0.0261.2 - 0.166 exp(-0.893r) rL, f 2 ( r ) = -0.0632 - 0.107r
(A2)
('43)
The Piare polynomials in the reduced Cartesian coordinates x = x'/L,, y = y'/L,, and z = .'Itx (x',y',z'the Cartesian coordinates): 4
P, = 1
+ I=c [ar(x2'+ y2') + bp2' + c,1(x21+ y2')z2'] 1
(A4)
4
P2 = 1
+ I=c [ d I ( x 2 '+ y2') + elz2'] 1
('45)
4
= 1
+ I=cI If,(x2' + y2') + g/z2q
P4 = 1
+ c [ h I ( x 2 '+ y2') + kIz2']
P3
(A61
4
I= 1
('47)
with the coefficients given in Table I. The fitting procedure was performed in two steps: (i) The potential energy I/ was fitted in the range of possible reduced coordinates (1x1, bl 5 0.5; Izl < 1.25) by using 714 points on a grid with 1x1, bl = 0, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2.0.3, 0.4, 0.5 and IzI = 0, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25. (ii) When satisfactory agreement between calculated and fitted pair interaction energies (eq A l ) was obtained, the analytical derivative of the potential function with respect to the three reduced Cartesian coordinates was calculated and the exact Ewald forces on the grid ( 3 X 714 points) were fitted to the resulting expression. This procedure yielded a standard deviation CT = 2.8 X relative to the force between two charges a t a distance of 1 A. Total electrostatic forces both from the exact Ewald summation and from the fitted relationship were calculated for the particles of three configurations of an equilibrated simulation of the same system (except for the use of the shifted force method by Streett et al.39for the Coulombic interactions). The average relative error for the fitted forces parallel and perpendicular to the surface was about 0.2% and 0.03%, respectively. These errors are smaller by an order of magnitude than the one given by he ye^^^ using q = L,/ 1.2 and evaluating the first eight replica cells only whereas the computational effort is still smaller. It should however be noted that the fit parameters depend on the elongation of the basic cell, especially because eq A1 does not have the correct asymptotic behavior for lzl m.
-
Registry No. H20, 7732-18-5; Pt, 7440-06-4.
