Zirconia Interface: Molecular

J. Phys. Chem. B 2004, 108, 19687-19697

19687

Solvation Dynamics at the Water/Zirconia Interface: Molecular Dynamics Simulations† Lucimara R. Martins and Munir S. Skaf* Instituto de Quı´mica, UniVersidade Estadual de Campinas, Cx. P. 6154, Campinas-SP 13083-970, Brazil

Branka M. Ladanyi* Department of Chemistry, Colorado State UniVersity, Fort Collins, Colorado 80523 ReceiVed: July 2, 2004

We report the results of a molecular dynamics (MD) simulation study of solvation dynamics associated with electronic excitation of the coumarin 343 (C343) dye at the interface of water and zirconia (ZrO2). We use an all-atom representation of the species present in this system. The ZrO2 partial charges and C343 geometries and charge distributions in the ground and excited electronic states have been obtained from electronic structure calculations. Our work is inspired by recent time-resolved fluorescence experiments involving C343 adsorbed at the surface of zirconia nanoparticles dissolved in H2O and D2O (Pant, D.; Levinger, N. E. J. Phys. Chem. B 1999, 103, 7846). In addition to simulation of solvation dynamics, we investigate the structure and dynamics of the dye and water in the presence of a planar ZrO2 interface. We also address several issues relevant to the interpretation of experiments, including the solvent isotope effects and the ionization state of the carboxylate group of C343 on solvation dynamics. We find that the neutral form (C343) of the dye is more strongly adsorbed at the ZrO2 interface and that the water portion of the solvation response for this form of the dye is significantly faster than for the deprotonated form, C343-. We also find that D2O-H2O solvent isotope effects on the solvation response of either form of the dye are quite modest and independent of the presence of zirconia. Rotational motion of the solute relative to the ZrO2 surface contributes a significant, very slowly relaxing component to the interfacial solvation dynamics. We discuss the implications of our findings for the interpretation of the experimental data.

I. Introduction Many technologically important processes occur at solid/ liquid interfaces. Understanding the structure and dynamics of liquids at solid interfaces impacts electrochemistry, heterogeneous catalysis, solar energy conversion, chromatography, fluid transport through clay minerals, and numerous other natural processes and technological applications. Time-resolved electronic spectroscopy using dyes that reside preferentially at or near the interface provides the means of probing solvation dynamics in heterogeneous environments. Several such studies involving solid/liquid,1-4 liquid/liquid,5,6 and liquid/vapor7-11 interfaces have been reported. In addition to experiments, molecular simulation designed to model solvation dynamics at several types of liquid interfaces have been carried out.12-17 While solvation dynamics in bulk polar liquids proceeds mainly through orientational relaxation of the surrounding solvent,18 additional solvation mechanisms can become important in heterogeneous systems.12-17,19-21 They include changes in the location of the probe molecule upon electronic excitation and structural rearrangements of the interface in the vicinity of the probe molecule. To what extent these additional mechanisms contribute is difficult to determine by purely experimental means, especially in view of the fact that interfacial dynamics of most liquids differs from its bulk counterpart. In this paper, we present molecular modeling and simulation results for solvation dynamics at a water/zirconia (ZrO2) †

Part of the special issue “Frank H. Stillinger Festschrift”. * Corresponding authors. E-mail addresses: [email protected]; [email protected].

interface. This study is motivated by the fact that understanding electron-transfer processes at such liquid/semiconductor interfaces is an important step in the design of photoelectrochemical cells22-24 and that adsorption of dyes and solvation dynamics accompanying their electronic excitation at such interfaces have not yet been studied by molecular simulation. We have chosen zirconia since a number of interesting experimental results on the properties of the coumarin 343 (C343) dye adsorbed on ZrO2 nanocrystals in aqueous solution are available2,3 and because it serves as a model for titania (TiO2), a promising medium for solar energy conversion.22-28 Pant and Levinger studied solvation dynamics of water at the zirconia (ZrO2) surface using time-resolved fluorescence spectroscopy.2,3 They found that interfacial solvent response displays two subpicosecond diffusive components with relaxation times similar to those of bulk aqueous solution. However, the relative amplitude corresponding to the shorter component was significantly larger for the ZrO2/water interface than for the bulk solution, which led to a faster average solvation response for molecules at the ZrO2 surface. Although a significant isotope effect is observed for the longer solvation time component in bulk D2O relative to H2O,29 no isotope effect was observed for the solvation dynamics at the surface of ZrO2.3 Surprisingly, all three characteristic time responses, (τ, τ1/e, and 〈τ〉), were found to be essentially the same for H2O and D2O at the ZrO2 surface. Several possible explanations were presented for the faster solvation dynamics and the missing isotope effect. First, the nanoparticle surface blocks a part of the dye molecule, so the relative number of solvent molecules surrounding the

10.1021/jp0470896 CCC: $27.50 © 2004 American Chemical Society Published on Web 09/25/2004

19688 J. Phys. Chem. B, Vol. 108, No. 51, 2004 C343 molecule is smaller than the number surrounding it in bulk solution. A smaller number of water molecules in the solvation shell of the dye could lead to a shorter relaxation time if collective solvent motion were responsible for the observed relaxation. Second, the water molecules at the surface could be prealigned. If the water molecules are close to the position required to solvate the excited state of the dye then we would anticipate a smaller degree of reorientation leading to a smaller Stokes shift. However, neither of these effects can account for the missing isotope effect. It was suggested that the most plausible reason for the faster response2 and the missing isotope effect3 at the surface is that the hydrogen bonding is disrupted allowing water to move more freely than it does in bulk solutions. Largely motivated by Pant and Levinger’s experimental results, we have undertaken molecular dynamics (MD) simulations of C343 at H2O/ZrO2 and D2O/ZrO2 interfaces to better understand at the molecular level the main features of solvation dynamics in these systems. The reported solvation dynamics experiments for bulk and interfacial systems have been carried out at different pH to solubilize the dye in the corresponding solutions.2,3 In the presence of the zirconia nanoparticles, the solutions are acidic, and the dye is found predominantly in its protonated form (C343). In contrast, the fluorescence measurements in bulk water are carried out for alkaline solutions, and the dye is predominantly deprotonated (C343-). With this in mind, we have performed simulations with both C343 and C343-. Simulations were carried out in both H2O and D2O to shed light on the isotope effects on interfacial and bulk solvation dynamics. Apart from the solvation responses themselves, which are directly comparable to the available experimental data, we have also investigated the structural and dynamical features of the solute and solvent molecules with respect to the zirconia interface. In section II, we report quantum chemical results for the parametrization of the ZrO2 crystal and the solute probes. Section III describes details of the simulations and interaction potential energies. Our results are presented and discussed in section IV. Section V contains a summary of our main findings and concluding remarks. II. Quantum Parametrizations A. Zirconia. Under standard conditions, the stable phase of pure zirconia (baddeleyite) forms a monoclinic crystal (m-ZrO2), in which each Zr is coordinated by seven oxygen atoms. With increasing temperature, the crystal structure transforms into tetragonal and then to cubic and, with increasing pressure, to orthorhombic polymorphs. The m-ZrO2 crystal structure has been well-characterized by neutron diffraction.30 The structure can be seen as alternating layers of Zr and O atoms where the Zr coordination varies between 3 and 4 from one layer of oxygens to the next. The oxygen atoms are divided into two main groups according to the Zr-O distances. The average ZrOI distance is 2.07 Å, while the average Zr-OII separation is 2.21 Å. The OII atoms form a planar square that corresponds to half of a cubic structure. The OI atoms form a triangle nearly parallel to the OII plane. The structure of the m-ZrO2 unit cell is shown in Figure 1. Monoclinic zirconia slabs of different dimensions were built using CRYSTAL9831 on which ab initio quantum chemical calculations were performed to obtain partial charges for the Zr and O atoms after proper addition of hydrogen atoms to the dangling bonds. Mulliken charges were obtained for a supercell containing eight unit cells (2 × 2 × 2) at the RHF/STO-3G

Martins et al.

Figure 1. ZrO2 unit cell and a (6 × 6 × 2) ZrO2 crystal slab used in the simulations inside the simulation box.

level. Since the charges on the OI and OII turned out very similar, we have made no distinction between the oxygen types and computed the mean O charge from the oxygens bound to the central Zr atoms only. The results are summarized in Table 1. B. Coumarins C343 and C343-. The parametrizations for both the protonated (C343) and anionic (C343-) dyes were performed according to the procedure used for coumarin 153, described in previous works.32,33 For the ground state, a full geometry optimization was carried out from ab initio calculations at the RHF/6-31G(d,p) level using Gaussian 98.34 The molecular geometry is also optimized in the excited state through a configuration interaction procedure with single excitations (CIS). For both ground and excited states, the atomic charges were derived from electrostatic potential fits using the MerzSingh-Kollman method.35-37 Due to the presence of the carboxylic group, there are several possible conformations in the electronic ground state. Our calculations show that the C343 structure depicted in Figure 2, with an intramolecular H-bond, represents the global minimum on the molecular energy surface, in close agreement with the results of Cave and Castner.38 These authors have also shown that only this conformation is thermally accessible in the ground state since the next stable conformer lies 5.84 kcal/mol higher in energy. The calculated ground-state dipole moment of C343 is 10.91 D, which compares well with the experimental estimates (9.86 D) by Moylan39 and calculations reported by Cave and Castner38 (11.76 and 10.66 D, using DFT and MP2, respectively). The dipole moments of the S1 and S0 states of C343 differ by 3.58 D, in close agreement with the TDDFT results of ref 38 (4.06 D). For C343-, the ground- and excited-state dipole moments at the present level of calculation are 24.21 and 25.17 D, respectively. The change in dipole moment (0.96 D) is less than half that reported by Song and Chandler40 using the semiempirical MNDO method. Such a difference stems from the orientation of the carboxylic group (COO-). In the MNDOoptimized structure, the COO- is perpendicular to the plane of the molecule, and the adjacent six-membered ring is planar. At the RHF level, the COO- lies closer to the conjugated plane with some distortions of the adjacent ring. The RHF energy

Solvation Dynamics at the Water/Zirconia Interface

J. Phys. Chem. B, Vol. 108, No. 51, 2004 19689

TABLE 1: Zirconia and Coumarin Lennard-Jones and Charge Parameters ZrO2 atom

σ (Å)

/kB (K)

q (e)

Zr O

5.714 2.630

0.014 140.136

0.7952 -0.3976

C343 (C343-)a atom

σ (nm)

/kB (K)

qS0 (e)

qS1 (e)

C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14 C15 H16 N17 O18 O19 H20 H21 H22 H23 H24 H25 H26 H27 H28 H29 H30 H31 H32 C33 O34 O35 H36

0.35 0.355 0.355 0.35 0.35 0.35 0.35 0.35 0.355 0.355 0.355 0.355 0.375 0.355 0.355 0.242 0.325 0.29 0.296 0.242 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.375 0.296 0.300 0.000

33.2 35.2 35.2 33.2 33.2 33.2 33.2 33.2 35.2 35.2 35.2 35.2 52.8 38.2 38.2 15.1 85.5 70.4 105.7 15.1 15.1 15.1 15.1 15.1 15.1 15.1 15.1 15.1 15.1 15.1 15.1 15.1 52.8 105.7 85.5 0.0

-0.067 (-0.073) -0.242 (-0.162) 0.208 (0.131) -0.116 (-0.078) -0.026 (-0.025) -0.112 (-0.059) -0.047 (-0.033) -0.057 (-0.029) -0.012 (-0.010) 0.498 (0.417) -0.280 (-0.135) -0.254 (-0.311) 1.008 (0.987) -0.720 (-0.452) 0.282 (0.013) 0.123 (0.146) -0.145 (-0.245) -0.477 (-0.532) -0.661 (-0.615) 0.198 (0.186) 0.088 (0.066) 0.084 (0.069) 0.042 (0.020) 0.042 (0.034) 0.076 (0.076) 0.069 (0.059) 0.093 (0.067) 0.079 (0.057) 0.048 (0.016) 0.041 (0.032) 0.051 (0.033) 0.052 (0.033) 0.931 (0.900) -0.601 (-0.813) -0.687 (-0.771) 0.491

-0.073 (-0.030) -0.214 (-0.289) 0.208 (0.131) -0.149 (-0.111) 0.009 (-0.017) -0.172 (-0.049) 0.0002 (-0.061) -0.123 (-0.028) 0.119 (0.055) 0.442 (0.549) -0.029 (-0.144) -0.397 (-0.386) 0.885 (0.872) -0.574 (-0.412) -0.021 (-0.029) 0.152 (0.159) -0.014 (-0.189) -0.452 (-0.533) -0.628 (-0.575) 0.217 (0.193) 0.094 (0.075) 0.093 (0.073) 0.039 (0.018) 0.037 (0.037) 0.080 (0.065) 0.072 (0.053) 0.106 (0.065) 0.092 (0.053) 0.049 (0.026) 0.036 (0.038) 0.070 (0.037) 0.070 (0.035) 0.896 (0.892) -0.593 (-0.798) -0.692 (-0.777) 0.478

a

Charges for C343- are in parentheses.

the numbering of the atoms relative to the molecular structure and the orientations and relative sizes of the ground- and excitedstate dipole moments of C343. III. Potentials and Simulation Details All molecules are treated as rigid bodies, represented by nonpolarizable interaction site models consisting of a sum of Lennard-Jones plus Coulombic terms between all pairs of sites. Lennard-Jones parameters for the zirconia atoms were obtained from the Buckingham potential determined by Minervini et al.41 according to the procedure described by Rappe´ et al.42 We have used σ ) 12 to achieve best agreement between the LennardJones and Buckingham potentials at long range, as suggested by Goddard and co-workers.43 The results are shown in Table 1. Lennard-Jones parameters from the OPLS-AA force field were used for the coumarin solutes, and the solvent is SPC/E water.44 Usual Lorentz-Berthelot combination rules are used for the Lennard-Jones interaction between different types of atoms.45 We have considered systems consisting of a single C343 or C343- solute, a (6 × 6 × 2) slab of zirconia crystal, and 642 SPC/E water molecules with cubic periodic boundary conditions. Initial configurations were obtained using the molecular packing program PACKMOL.46 The system size was chosen after preliminary tests to ensure that boundary effects are negligible, especially in the Z direction (normal to the water/zirconia interface). The simulations were performed in the NVE ensemble at an average temperature of 298 K. Lennard-Jones forces were cut off at half the box length, and long-range Coulombic interactions were treated using Ewald sums. A neutralizing background was added in systems containing the C343- solute. The equations of motion for the solute and water molecules were solved using SHAKE and leapfrog algorithms45 with a time step of 3 fs. The relative positions of Zr and O atoms in zirconia were kept fixed at the crystal structure throughout the simulations. Long (∼3 ns) equilibrium NVE production runs were generated with occasional intervals of velocity rescaling to maintain the desired temperature. No trajectory data were collected during rescaling. The thermalization procedure and equilibrium production runs were performed for the protonated and anionic dyes for both S0 and S1 state charge distributions in the water/zirconia system and in neat bulk water. Similar runs were also performed for interfacial and bulk D2O systems. Nonequilibrium solvation dynamics runs were generated in the usual way by choosing a series of statistically uncorrelated configurations from the equilibrium runs with the solute in the S0 state as starting points of nonequilibrium trajectories. At t ) 0, the solute charge distribution was switched to that corresponding to the S1 state, and the system was allowed to relax under the new solute-solvent interactions. To generate reasonable statistics for the solvation response, we had to perform a fairly large number of such runs (400), each lasting 15 ps of real dynamics. IV. Results and Discussion

Figure 2. C343 molecular structure and atom numbering. The arrows represent the directions and relative magnitudes of the molecular dipoles in the S0 (s) and S1 (- - -) states.

difference between the planar and perpendicular -COOorientations is 24 kcal/mol. The partial charges for the S0 and S1 states of protonated and deprotonated coumarins are listed in Table 1. Figure 2 shows

A. Coumarin Adsorption. It has been suggested that the C343 location is one of the reasons for the observed faster solvation dynamics at the interface in comparison to the bulk solution. If the C343 molecule is adsorbed onto the nanoparticle surface, a portion of it will be inaccessible to the solvent molecules, so the number of water molecules surrounding the C343 solute will be smaller than in bulk solution.3 To analyze this possibility, we have calculated the distribution of distances between the solute (either C343 or C343-) center

19690 J. Phys. Chem. B, Vol. 108, No. 51, 2004

Figure 3. Locations of C343 and C343- relative to the zirconia surface: (a) probability distribution for the coumarin center-of-mass distance from the zirconia surface; (b-d) angle probability distributions for different coumarin body-fixed vectors relative to the surface normal. Bottom drawings show the average orientations for the ground and excited states of C343 adsorbed onto the zirconia surface.

of mass and the zirconia interface. The probability distribution functions are shown in Figure 3. The protonated C343 solute lies at an average distance of 2.8 Å from the crystal surface for both the S0 and S1 states and remains in the vicinity of the interface during the entire course of the simulation. In the S1 state, the distribution is a little broader on the right side of the peak, corresponding to larger Z values, indicating that the center of mass of the solute in the excited state may be found at somewhat larger distances. For comparison, the distance distributions for the deprotonated dye in S0 and S1 states are also shown. The peaks of the C343- distributions are located a little further away from the interface and are broadened toward larger Z values, indicating that C343- may wander further away from the zirconia surface. Experiments suggest that C343- is adsorbed more weakly than C343, in agreement with this finding. However, like the protonated form, C343- in our model seems to be predominantly found in the interfacial region. It is possible that including anion polarizability into the model for C343would increase its tendency to migrate into bulk water. The relative orientation of the solute with respect to the interfacial plane is also important since it can influence how exposed the dye is to the solvent molecules. The orientation of the solute relative to the interface normal (Z direction) has been determined by computing the probability distributions for the angles θpl, θperp, and θdip, between Z and three characteristic unit vectors: parallel (upl) and perpendicular (uperp) to the molecular (fused aromatic ring) plane and along the dipole moment (udip). More specifically, upl is along the vector linking the atoms C2 and C12, while udip is close to being in the molecular plane and in the case of C343, approximately perpendicular to upl (cf. Figure 2). With peaks at θpl = 90° and

Martins et al. θperp = 10°, the angle distributions shown in Figure 3 indicate that the C343 molecule is almost parallel to the zirconia surface. This means that ground state of C343 is surrounded by water essentially only on one side of the molecular group plane. Water access to the other side is blocked by the interface. In the excited state, C343 has a significantly different average orientation. For S1, the peak of the distribution is at θpl = 120°, indicating that the molecular plane is not parallel to interface. For both θpl and θperp, the angle distributions are also broader, suggesting that the solute in the excited state may swing more widely relative to the interface but keeping the dipole vector oriented approximately in the same direction, given that the distribution of θdip remains almost the same for C343 in the two electronic states. The drawings at the bottom of Figure 3 depict the C343 molecule in its average position relative to the interface for the S0 and S1 states. At the average S1 orientation, the C6, C7, C8, C9, C11, C12, and C15 sites (see Figure 2 for the site labels) have moved further away from the interface compared to where they were in the S0 state, whereas the carboxylic groups remain close to the zirconia surface. The average number of solvent molecules contained within the nonspherical first coordination shell of the solute in the S0 and S1 states has been computed by measuring distances between the solvent center of mass and the nearest solute atom and using the concept of “solvation shell distribution” introduced by Song et al.47 in their study of supercritical solvation of molecules of arbitrary shape. Detailed definition of the shell distribution and how to compute the coordination number within the first nonspherical solvation shell are given in ref 47. From the equilibrium simulations, we obtain 38.4 and 38.6 water molecules within a shell 5 Å thick surrounding the solute for the S0 and S1 states, respectively. This indicates that the coordination number is independent of the electronic state of the solute within our parametrization of the dye. However, using the same shell thickness for C343 immersed in bulk water, we find a coordination number of 53 water molecules, which is consistent with the discussion above. B. Coumarin Dynamics. The translational dynamics of the coumarin in the presence of the interface is investigated by examining of the mean-squared displacements, 〈[∆r(t)]2〉, computed separately for the X, Y, and Z directions as shown in Figure 4. As a consequence of C343 adsorption onto the zirconia surface, the total mean-squared displacement in the 10 ps time interval is roughly 5 times smaller than in the bulk. The presence of the solid surface also affects differently the displacements in the perpendicular (Z) and parallel (XY) directions, the former being much smaller on the time scale of our observations, indicating that C343 translates mainly in the plane of the ZrO2 interface. Figure 4 also shows that translational motions in the bulk and at the interface are overall less hindered for the protonated dye due to the stronger electrostatic interactions between water and C343-. Note, however, that the anionic form of the dye translates somewhat more easily in the direction perpendicular to the interface than does C343. This is consistent with the fact that C343- is more weakly adsorbed, as shown in Figure 3a. The time correlation functions (TCFs),

C2(t) ) 〈P2[u(0)‚u(t)]〉

(1)

where P2(x) ) (3x2 - 1)/2 is the second-order Legendre polynomial, for the upl, uperp, and udip vectors are shown in Figure 5 for the chromophore in bulk water and at the ZrO2/ H2O interface. Given that the orientations of udip in the S0 and


J. Phys. Chem. B, Vol. 108, No. 51, 2004 19691 C. Solvation Structures. The equilibrium solvation structures for the ground- and excited-state C343 are not very different. Thus, small solvent structural rearrangements are involved in the solvation dynamics. There are only small density variations in the vicinity of some solute sites after excitation of the probe. Overall, a slight increase is found in the solvent density around the atoms C6, C7, C8, C9, C11, C12, and C15, located at the bottom of the molecule (cf. Figure 2) and a small decrease around the upper part of C343 (sites C1, C5, C13, O18, and O19). These small rearrangements are associated with the changes in the relative orientation of C343 with respect to the interface upon excitation, but not so much with the charge jumps on each site individually. The solvation structures near sites C2, C3, C4, C10, C14, and N17 are almost invariant with the excitation of the solute. In Figure 6, the radial distribution functions between water and the sites C15, C13, and N17 are shown as examples of increased, decreased, and invariant solvation structures, respectively. Nonequilibrium solvation structures were also computed from the nonequilibrium runs. However, no noticeable changes in the radial distribution functions in comparison to the S0 structures were observed within the time frame of these runs (∼15 ps). D. Solvation Dynamics. We now focus on the observables that are directly accessible from dynamic fluorescence Stokes shift measurements. The key quantity is the solvation response function, S(t), which is obtained from time-resolved fluorescence spectroscopy by monitoring the evolution of the position of the emission band according to

Figure 4. Mean-squared displacements for C343 (s) and C343(- - -) in bulk (top panel) and interfacial (bottom panel) systems.

TABLE 2: Solute P2 Orientational Relaxation Times in Bulk H2O and at the H2O/ZrO2 Interface Obtained from TCF Slopes in the 0.4 < t < 12.0 ps Time Interval system

τpl (ps)

τperp (ps)

τdip (ps)

C343 in bulk C343- in bulk C343 at interface C343- at interface

31.1 39.5 224 244

32.3 39.2 403 345

46.6 179 505 2877

S1 states of C343 are almost the same, the characteristic time associated with the decay of C2,dip(t) is likely to be closely related to the fluorescence depolarization decay time, which has been measured experimentally.2,3 The decay times for all three vectors, estimated by assuming that the decay is exponential after about 1 ps, are summarized in Table 2. It can be seen from Figure 5 that the decay rate of C2(t) in the interfacial system is much slower than in the bulk. It is also more anisotropic: In the bulk, the decay rates of C2,pl(t) and C2,perp(t) are essentially the same, while at the interface C2,perp(t) decays considerably more slowly than does C2,pl(t). This is consistent with the fact that the molecule is adsorbed with the aromatic ring plane parallel to the zirconia surface and indicates that reorientation in this plane is less hindered than out-of-plane rotational relaxation. When orientational relaxation rates of C343 and C343- are compared, C2,dip(t) is more slowly relaxing for the anionic form of the solute, reflecting stronger effects of dielectric friction. The characteristic dipole reorientational times of C343 and C343- in the interfacial system, estimated the slope of ln C2,dip(t), are about 0.5 and 2.9 ns, more than an order of magnitude longer compared with their bulk counterparts of 47 and 180 ps but somewhat lower than 5.7 ns measured experimentally.3 The interfacial MD characteristic times have to be interpreted with caution since they are much longer than the time interval over which the time-correlation functions were calculated.

S(t) )

hν(t) - hν(∞) hν(0) - hν(∞)

(2)

where ν(t) is the maximum of the fluorescence band at time t after the probe excitation at t ) 0 and ν(∞) corresponds to the steady-state emission frequency. The solvation response function can be calculated from the nonequilibrium MD trajectories within the Franck-Condon approximation by computing the difference between solute-solvent electrostatic potential energies, USel1 and USel0, for the solute in electronically excited (S1) and ground (S0) states:

∆E(t) ) USel1(t) - USel0(t)

(3)

In terms of ∆E, S(t) can be expressed as

S(t) )

〈∆E(t)〉ne - 〈∆E〉S1 〈∆E〉S0 - 〈∆E〉S1

(4)

where 〈...〉ne indicates an average over nonequilibrium trajectories, while 〈...〉S0 and 〈...〉S1 denote equilibrium averages for systems that contain, respectively, the ground (S0) and the excited (S1) state solute. Since ∆E is pairwise-additive in our model, the separate water and zirconia contributions to the total response can be calculated from the water-coumarin and zirconia-coumarin interactions:

S(t) ) SH2O(t) + SZrO2(t)

(5)

where the response for the coumarin interaction with solvent component R (H2O, ZrO2)

SR(t) )

〈∆ER(t)〉ne - 〈∆ER〉S1 〈∆E〉S0 - 〈∆E〉S1

(6)


Martins et al.

Figure 5. Reorientational TCFs, C2(t), for C343 (left panels) and C343- (right panels) dyes in bulk (top panels) and interfacial (bottom panels) environments. Different line styles correspond to TCFs for different unit vectors, upl, uperp, and udip, embedded in the solute molecules.

Results for the total response and its separate contributions are shown in the upper panel of Figure 7 for the protonated dye. SH2O(t) is responsible for roughly 90% of the initial response and is completely relaxed within 4-6 ps. The remaining 10% of the initial response is due to the coupling between C343 and ZrO2, which relaxes on a much longer time scale. Given that the zirconia slab is stationary in our model, the relaxation of SZrO2(t) can occur only through motion of C343 relative to the interface. Thus the slow reorientational dynamics of C343 adsorbed on the zirconia surface leads to a slowly relaxing component in the total S(t). Despite the zirconia contribution being small, it determines entirely the solvation response for t > 5 ps. Even 15 ps after excitation, the total response is still roughly 8% from being completely relaxed. In the lower panel of Figure 7, we compare the interfacial SH2O(t) with the solvation response for C343 in bulk water. The water response in the interfacial system is slightly slower than the response of bulk water (the differences between SH2O(t) and Sbulk(t) seen in Figure 7 would increase a little if SH2O(t) were normalized to 1 at t ) 0). In Figure 8, we compare the interfacial S(t) for C343 to its linear-response counterparts, the solvation TCFs,48

C0,1(t) )

〈δ∆E(t)δ∆E(0)〉S0,S1 〈(δ∆E)2〉S0,S1

(7)

where δ∆E(t) ) ∆E(t) - 〈∆E(0)〉S0,S1 is the fluctuation of the solute-solvent energy gap computed from the equilibrium runs in the S0 or S1 states. In highly polar solvents, including water,49-52 S(t) is often found to agree well with the solvation TCFs.

As can be seen from Figure 8, this is true in the present case: C0(t) and C1(t) are nearly identical, and both are in good agreement with S(t), indicating that the linear response approximation is reasonable in this case. In view of this and given that high-quality data for the solvation TCFs can be obtained with less computational effort than for S(t), we use C0(t) instead of the noisier S(t) in the subsequent analysis of the solvation response. In Figure 9, we display the comparison with experimental data for solvation dynamics in bulk water (upper panel) and at the H2O/ZrO2 interface. In the case of the interfacial system, the experimental response was observed to essentially level off at 5 ps, and the remainder of relaxation, which occurs essentially on the time scale of solute reorientation (around 5.7 ns), was not detected in the time scans extending to about 500 ps.3 In view of this fact, we compare the experimental solvation response to the normalized [C0(t) - C0(5ps)] in the case of this system, while the full C0(t) is used in the case of solvation dynamics in bulk water. The data shown in Figure 9 indicate that the MD results are in very good agreement with experiment. The short-time librational oscillations that are observed in the calculated C0(t)’s are below the experimental time resolution, so they cannot be expected to appear in the experimental S(t)’s. Note that the MD data are for C343- in the bulk and for C343 at the interface. In Table 3, we summarize the results of multiexponential fits to our C0(t) data and also include the experimental S(t) data in the same format. In addition to solvation dynamics results, we have included the results for the steady-state Stokes shifts,

∆ν ) [〈∆E〉S0 - 〈∆E〉S1]/(hc)

(8)


Figure 6. Equilibrium site-site solute-water pair distribution functions in the presence of zirconia interface. Solid and dashed lines are for the ground and excited states of C343, respectively.

The experimental ∆ν values are larger than the calculated ones, suggesting that the change in the solute-solvent electrostatic interactions associated with the S0 f S1 transition is underestimated in our calculations. It is likely that interactions with water lead to polarization of the solute and an increase in its effective dipole, as has been observed in the case of solvation of coumarin dyes in other polar solvents.53 We expect that addition of such interactions would improve agreement with experimental ∆ν data in the present system. Figure 10 illustrates the effects of the ionization state of the dye on solvation dynamics. Again, the top panel depicts the results for the bulk solvent (H2O), while the bottom panel is for the H2O/ZrO2 interface. The results show that solvation dynamics of the protonated C343 solute is faster than that of the deprotonated species, both in bulk and at the interface. In fact, the solvation dynamics of C343 in H2O/ZrO2 is faster than the solvation dynamics of C343- in bulk H2O for t < 5 ps (for longer times, the zirconia contribution starts to slow C0(t)). This is a very interesting result given that in the solutions with the zirconia nanoparticles where the experiments have been carried out, the dye is predominantly protonated.2,3 On the other hand, in bulk solutions the C343- anionic species prodominates. Therefore, simulations suggest that the ionization state of the dye may be a key factor underlying the observed faster solvation dynamics in the presence of the zirconia interface. Figure 11 illustrates the MD results for the solvent isotope effects on the solvation response. The upper panel depicts C0(t) for C343- in bulk H2O and D2O, while the lower panel illustrates the behavior of C0(t) for C343 in H2O/ZrO2 and in D2O/ZrO2. It can be seen that solvation dynamics slows slightly when D2O is substituted for H2O and that the effect of this


Figure 7. Nonequilibrium solvation responses of C343 at the H2O/ ZrO2 interface and in bulk water: (top panel) components of the response at the interface; (bottom panel) comparison of pure and interfacial water solvation responses. The insert shows data for short times.

Figure 8. Comparison of the nonequilibrium solvation response, S(t), and the solvation TCFs, C0(t) and C1(t), for the C343 solute in its S0 and S1 states at the H2O/ZrO2 interface.

substitution is about the same in the bulk liquid and at the interface. Thus our simulations indicate that if a single form of the solute (C343- or C343) is present, there is very little difference between the isotope effects in the bulk and at the interface. It is worth noting in this regard that there is a significant H2O/D2O isotope solvent effect on the pKa’s of carboxylic acids.54 The dissociation constant in D2O is several times smaller than that in H2O, which would result in a larger relative concentration of the neutral form of the dye in the deuterated solvent. The fact that solvation dynamics is faster for C343 than


Martins et al.

Figure 9. Comparison with experiment for solvation dynamics of C343- in bulk water (top panel) and C343 at the H2O/ZrO2 interface (bottom panel). In the case of solvation in bulk liquid, the experimental response is compared to C0(t). In the case of solvation in H2O/ZrO2, the experimental response is compared to [C0(t) - C0(5ps)]/[1 - C0(5ps)].

TABLE 3: Multiexponential Fit Parameters to the Linear Response Solvation Dynamics, C0(t), from the S0 Equilibrium Runs and Stokes Shifts (in cm-1)a a1

τ1

a3

τ2

a3

Experimental C343- in H2O 0.34 0.07 0.66 0.63 C343- in D2O 0.18 0.06 0.82 1.01 C343 in 0.73 0.09 0.27 0.64 H2O/ZrO2 C343 in 0.63 0.06 0.37 0.56 D2O/ZrO2 C343 in H2O C343 in D2O C343- in H2O C343- in D2O C343 in H2O/ZrO2 C343 in D2O/ZrO2 C343- in H2O/ZrO2 C343- in D2O/ZrO2 a

0.76 0.76 0.63 0.64 0.61

0.06 0.08 0.09 0.10 0.05

0.24 0.24 0.37 0.36 0.26

MD 0.78 1.09 1.14 1.45 0.45 0.13

τ3

τ0

〈τ〉

τ1/e

∆ν

0.17 0.44 0.40 800 0.26 0.84 0.81 500 0.12 0.24 0.14 150 0.09 0.25 0.12 190

0.08 0.10 0.14 0.15 5.6 0.08

0.23 0.32 0.48 0.59 0.88

0.07 0.10 0.22 0.23 0.13

414 427 132 131 394

0.65 0.06 0.25 0.73 0.10 12.6 0.09 1.48 0.14 402 0.56 0.08 0.37 1.42 0.07 13.8 0.14 1.53 0.35 125 0.58 0.09 0.31 1.41 0.11 13.1 0.15 1.93 0.34 126

Times are given in ps. Experimental data are from ref 3.

for C343- would compensate to some extent for the slower rate of D2O solvent relaxation and might explain the diminished isotope effect on solvation dynamics at the zirconia interface. Thus the displacement by the solvent of the ionization equilibrium appears to be the likely candidate for the observed experimental trends within our model of solute-solvent interactions.

Figure 10. The effect of the solute ionization state on solvation dynamics in bulk water (top panel) and at the H2O/ZrO2 interface (bottom panel). Depicted are solvation TCFs, C0(t), for C343 (s) and C343- (- - -).

E. Structure and Dynamics of Water in the Presence of the ZrO2 Interface. In this section, we characterize the solvent structure and dynamics near the interface. To proceed with the analyses, we have defined distinct layers from the interface according to the distribution P(Z) of water oxygen-interface distances. Distribution functions for the H2O/ZrO2 system in the absence (solid line) and in the presence (dashed line) of the C343 solute are shown in Figure 12. The distributions are similar, but there are noticeable differences due to the solute, especially for Z > 5 Å. The overall trend, which should disappear for systems containing a smaller solute/solvent mole ratio, is that water is displaced by C343 from the interfacial region and therefore has a higher probability of being located at larger Z. P(Z) exhibits a fairly high first peak and a smaller second peak, indicating that there is some attraction between water and zirconia. The first two local minima occur at 3.5 and 7.0 Å and delimit, respectively, the first and second solvent layers parallel to the zirconia crystal surface. The slab between Z ) 3.5 and 7.0 Å contains molecules that participate directly in the solvation of the dye, as well as molecules that are far from the solute, in the XY plane. Water molecules beyond 7 Å from the interface are considered bulk water. Figure 13 shows the normalized probability distribution for the angle between the water dipole and the interface normal for molecules in the bulk region (Z > 7 Å, solid line) and in the first and second interfacial layers with Z < 3.5 and 3.5 < Z < 7.0 Å, respectively. For an isotropic environment such as bulk liquids, the orientation of any molecular vector follows a uniform distribution, which in polar coordinates is given by P(θ) ) 1/2 sin θ. Indeed, the angle distribution in the bulk is sinusoidal. The distribution for molecules in the second layer is very close to the truly isotropic behavior. Small differences appear due to the interactions with the solute and neighboring



Figure 13. The probability distribution for the angle between water dipoles and the interface normal for several regions, as indicated.

TABLE 4: H-Bonding Statistics for Different Regions from the Interfacea region or system

f1

f2

f3

f4

f5

Z < 3.5 Å 3.5 < Z < 7 Å bulk (Z > 7 Å) pure water

0.040 0.008 0.007 0.008

0.244 0.075 0.075 0.076

0.451 0.320 0.319 0.315

0.239 0.517 0.517 0.527

0.024 0.078 0.079 0.073

a

Figure 11. The solvent isotope effect on solvation dynamics in bulk liquid (top panel) and at the zirconia interface (bottom panel). Depicted are solvation TCFs, C0(t), for C343- in bulk H2O and D2O (top panel) and for C343 at H2O/ZrO2 and D2O/ZrO2 interfaces (bottom panel).

Figure 12. The probability distribution, P(Z), for the distance between water molecules and the zirconia surface.

water molecules. Near the interface, however, there is preferential alignment of the dipoles parallel to the zirconia plane, typical of hydrophobic hydration.55 It is also consistent with the expectation that water hydrogens will be closer to the zirconia surface, given that the negatively charged oxygen sites of ZrO2 are more accessible (cf. Figure 1). Figure 13 also shows the angle distribution for the first layer of a water/zirconia interface in the absence of the solute. Another important structural feature is the statistics of hydrogen bonding in the different interfacial regions, which have been computed using a geometric criterion56,57 and are shown in Table 4. Water molecules located in the second layer and bulk regions exhibit a H-bond distribution practically identical

The quantities fi refer to the fraction of molecules making i H-bonds.

to that of pure water. From bulk to the interface region, there is an increase in the fraction of water molecules participating in one, two, or three H-bonds and a decrease in the fraction of molecules engaged in four or more H-bonds. The average number of H-bonds per molecule is reduced from 3.6 in the bulk to approximately 3.0 close to the interface. Such a reduction is not unexpected given the presence of zirconia. It has been suggested that one of the reasons for the missing isotope effect at the interface and faster interfacial solvation dynamics is that the hydrogen-bonding is substantially disrupted near the interface, allowing water to move more freely than it does in bulk solutions.3 Although the typical tetrahedral hydrogen-bonding network of bulk water is disrupted at the interface, the simulations indicate that interfacial water does not move faster than bulk water. The mean-squared displacements of the solvent in different interfacial layers have been computed for water/zirconia systems with and without the coumarin for both the ground and excited states and for both the ionization states of C343. The results are practically indistinguishable for all these solute forms, and hence, data are shown in Figure 14 only for protonated C343 in the S0 state. The mean-squared displacements of water molecules in different layers (Figure 14) show reduced translational mobility for interfacial water. Further away from the interface, the slope of the mean-squared displacement approaches that of bulk water. The slowing of water translational diffusion near the interface occurs not only in the Z direction, perpendicular to the surface but also in the XY plane. The reduction in translational mobility indicates that there is significant electrostatic attraction between water and the zirconia surface. In the case of purely hydrophobic hydration, the translational mobility of water near the interface typically increases.55,58 In addition, we have also investigated the translational dynamics of D2O. One can observe from Figure 14 that there is an isotope effect on mean-squared displacement for all regions considered. The slopes of these lines divided by 6, which correspond approximately to the local translational diffusion


Martins et al. results are also summarized in Table 5. The single-molecule dipole orientational relaxation time in the region Z > 7 Å is similar to its value in bulk SPC/E water61 and increases considerably as one approaches the interface. The H2O/D2O ratio of relaxation times is again smaller than the experimental ratio in the bulk65 but practically the same for bulk and interfacial water. In summary, the isotope effects in the simulated systems are essentially independent of the location of the water molecules relative to the interface. V. Concluding Remarks

Figure 14. Mean-squared displacements, 〈(∆r)2〉, for H2O and D2O in three interfacial layers. In the case of H2O, mean-squared displacements in directions parallel (X or Y) and perpendicular (Z) to the interface are also shown.

TABLE 5: Water Translational and Rotational Mobilities: Slopes (divided by 6) of the Solvent Mean-Squared Displacementsa and Single-Molecule Dipole Reorientational Timesb for Different Regions from the Interface d〈(∆r)2〉/dt (10-5 cm2/s)

1/ 6

region or system H2O D2O H2O/D2O Z < 3.5 Å 3.5 < Z < 7 Å bulk (Z > 7 Å) pure water

1.04 2.01 2.35 2.53

0.90 1.73 2.02 2.15

1.16 1.16 1.16 1.18

single-molecule dipole reorientation times (ps) H2O

D2O

10.44 12.31 5.19 6.26 4.66 5.52 4.38 5.27

D2O/H2O 1.18 1.21 1.19 1.20

a Computed over a 0.5 < t < 10.0 ps time interval. b Defined as time integrals of 〈u(0)‚u(t)〉. Fits to exponential decay in the time interval 0.2 < t < 15.0 ps are used to extrapolate long-time behavior.

tensor components,59,60 are listed in Table 5, along with the simulated and experimental diffusion coefficients of neat H2O and D2O. The ratio of diffusion coefficients or diffusion tensor components for these species indicates that the MD isotope effect on the translational dynamics of interfacial water is very similar to that of bulk water, which had previously been determined for SPC/E water isotopes by Svishchev and Kusalik.61 However, the simulations underestimate the isotope effect in real water, which is not surprising given that quantum effects associated with the isotopic substitution,62-64 such as the stiffening of the H-bond network, are ignored in these simulations. Similarly, we have also investigated the dipole reorientational dynamics for H2O and D2O in different interfacial layers. These

We have presented results of MD simulations of solvation dynamics accompanying electronic excitation of the C343 dye at the water/zirconia interface. We have also studied the translational and rotational mobilities of the dye in the presence of the interface as well as its average location relative to the zirconia surface. Interfacial properties of water were studied as well. The D2O-H2O isotope effects on solvation dynamics in bulk and interfacial water were determined. Since both the neutral, protonated form of the dye, C343, and the deprotonated form, C343-, are present in aqueous media, the electronic structures and gas-phase equilibrium geometries were determined for both forms and used in our MD calculations. We summarize here our main findings. Both C343 and C343- in their ground electronic states adsorb onto zirconia, with the molecular plane parallel to the crystal surface. In the S1 state, the dye partially desorbs by rotating into the aqueous phase. The protonated form of the dye adsorbs more strongly than does C343-. Solvation dynamics at the interface proceeds via two distinct mechanisms occurring on very different time scales. About 90% of the response can be attributed to the reorganization of water and occurs on a fast time scale very similar to that of bulk aqueous solvation. The remainder of S(t) relaxes through solute reorientation relative to the zirconia surface and occurs on a nanosecond scale. The total solvation response is predicted accurately using the linear response approximation. There is a larger difference in the solvation responses of C343 and C343-, the latter being significantly slower both in bulk and in interfacial systems. We find that the solvation dynamics of C343 in H2O/ZrO2 is faster than that of C343- in bulk H2O, thus suggesting that the ionization state of the dye may be a key factor in the experimentally observed faster solvation dynamics in the presence of zirconia. The solvent isotope effects accessible via classical MD simulation are relatively modest and about the same at the interface and in bulk solvent for a given form of the dye. Thus our tentative explanation of the smaller isotope effect at the interface than in the bulk that is observed experimentally is that a larger proportion of the dye is in its neutral form in the presence of D2O, given that the carboxylic acid dissociation constant is several times smaller in D2O than in H2O.54 The translational and rotational mobilities of water are reduced near the interface, an indication of attractive waterzirconia interactions. However, this reduction in mobility does not affect significantly the portion of the solvation response due to water. Our solvation dynamics results are in good agreement with experiment and our findings on the reduction in the solute rotational mobility in the presence of the zirconia interface are also consistent with experimental results. However, the steadystate Stokes shifts in the solute electronic spectra are underestimated in our model. Solute polarization by the surrounding medium and modification of the solute electronic structure due

Solvation Dynamics at the Water/Zirconia Interface to hydrogen bonding to water molecules are the likely sources of this discrepancy. We are in the process of carrying out quantum mechanical calculations that will help us develop models that incorporate these effects. Acknowledgment. It is a pleasure to dedicate this article to Frank Stillinger in honor of his seminal contributions to the statistical mechanics of condensed phase systems. The authors are grateful to Nancy E. Levinger for many stimulating discussions. L.R.M. thanks the Brazilian Agency CAPES for supporting her visit to CSU and FAPESP (Grant No. 00/048796) for a graduate fellowship. This work has been supported by grants from FAPESP (No. 01/09374-3 to M.S.S.) and from the U.S. National Science Foundation (No. CHE 9981539 to B.M.L.). References and Notes (1) Yanagimachi, M.; Tamai, N.; Masuhara, H. Chem. Phys. Lett. 1992, 200, 469. (2) Pant, D.; Levinger, N. E. Chem. Phys. Lett. 1998, 292, 200. (3) Pant, D.; Levinger, N. E. J. Phys. Chem. B 1999, 103, 7846. (4) Shang, X. M.; Benderskii, A. V.; Eisenthal, K. B. J. Phys. Chem. B 2001, 105, 11578. (5) Bessho, K.; Uchida, T.; Yamauchi, A.; Shioya, T.; Teramae, N. Chem. Phys. Lett. 1997, 264, 381. (6) Benderskii, A. V.; Eisenthal, K. B. J. Phys. Chem. A 2002, 106, 7482. (7) Zimdars, D.; Dadap, J. I.; Eisenthal, K. B.; Heinz, T. F. Chem. Phys. Lett. 1999, 301, 112. (8) Zimdars, D.; Eisenthal, K. B. J. Phys. Chem. A 1999, 103, 10567. (9) Zimdars, D.; Eisenthal, K. B. J. Phys. Chem. B 2001, 105, 3993. (10) Benderskii, A. V.; Eisenthal, K. B. J. Phys. Chem. B 2000, 104, 11723. (11) Benderskii, A. V.; Eisenthal, K. B. J. Phys. Chem. B 2001, 105, 6698. (12) Benjamin, I. J. Chem. Phys. 1991, 95, 3698. (13) Benjamin, I. Chem. Phys. 1994, 180, 287. (14) Michael, D.; Benjamin, I. J. Chem. Phys. 2001, 114, 2817. (15) Rose, D. A.; Benjamin, I. J. Phys. Chem. 1992, 96, 9561. (16) Vieceli, J.; Benjamin, I. J. Phys. Chem. B 2003, 107, 4801. (17) Chandra, A.; Senapati, S.; Sudha, D. J. Chem. Phys. 1998, 109, 10439. (18) Maroncelli, M. J. Mol. Liq. 1993, 57, 1. (19) Faeder, J.; Ladanyi, B. M. J. Phys. Chem. B 2001, 105, 11148. (20) Pal, S.; Balasubramanian, S.; Bagchi, B. J. Chem. Phys. 2004, 120, 1912. (21) Thompson, W. H. J. Chem. Phys. 2004, 120, 8125. (22) Kamat, P. V. Chem. ReV. 1993, 93, 267. (23) Hagfeldt, A.; Gratzel, M. Acc. Chem. Res. 2000, 33, 269. (24) Cahen, D.; Hodes, G.; Gra¨tzel, M.; Guillemoles, J. F.; Riess, I. J. Phys. Chem. B 2000, 104, 2053. (25) Rehm, J. M.; McLendon, G. L.; Nagasawa, Y.; Yoshihara, K.; Moser, J.; Gra¨tzel, M. J. Phys. Chem. 1996, 100, 9577. (26) Cherepy, N. J.; Smestad, G. P.; Gra¨tzel, M.; Zhang, J. Z. J. Phys. Chem. B 1997, 101, 9342. (27) Bauer, C.; Boschloo, G.; Mukhtar, E.; Hagfeldt, A. J. Phys. Chem. B 2002, 106, 12693. (28) Boschloo, G.; Hagfeldt, A. Chem. Phys. Lett. 2003, 370, 381. (29) Barbara, P. F.; Walker, G. C.; Kang, T. J.; Jarzeba, W. Proc. SPIEInt. Soc. Opt. Eng. 1990, 1209, 18. (30) Winterer, M.; Delaplane, R.; McGreevy, R. J. Appl. Crystallogr. 2002, 35, 434.

J. Phys. Chem. B, Vol. 108, No. 51, 2004 19697 (31) Saunders, V. R.; Dovesi, R.; Roetti, C.; Causa, M.; Harrison, N. M.; Orlando, R.; Zicovich-Wilson, C. M. CRYSTAL98, 1.0 ed.; University of Turin: Turin, Italy, 1998. (32) Martins, L. R.; Tamashiro, A.; Laria, D.; Skaf, M. S. J. Chem. Phys. 2003, 118, 5955. (33) Martins, L. R.; Skaf, M. S. Chem. Phys. Lett. 2003, 370, 683. (34) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision A.7; Gaussian, Inc.: Pittsburgh, PA, 1998. (35) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129. (36) Besler, B. H.; Merz, K. M.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431. (37) Cieplak, P.; Kollman, P. J. Comput. Chem. 1991, 12, 1232. (38) Cave, R. J.; Castner, E. W. J. Phys. Chem. A 2002, 106, 12117. (39) Moylan, C. R. J. Phys. Chem. 1994, 98, 13513. (40) Song, X.; Chandler, D. J. Chem. Phys. 1998, 108, 2594. (41) Minervini, L.; Grimes, R. W.; Sickafus, K. E. J. Am. Ceram. Soc. 2000, 83, 1873. (42) Rappe´, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (43) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897. (44) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (45) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford: New York, 1987. (46) Martinez, J. M.; Martinez, L. J. Comput. Chem. 2003, 24, 819. (47) Song, W.; Biswas, R.; Maroncelli, M. J. Phys. Chem. A 2000, 104, 6924. (48) Carter, E. A.; Hynes, J. T. J. Chem. Phys. 1991, 94, 5961. (49) Maroncelli, M.; Fleming, G. R. J. Chem. Phys. 1988, 89, 5044. (50) Bader, J. S.; Chandler, D. Chem. Phys. Lett. 1989, 157, 501. (51) Levy, R. M.; Kitchen, D. B.; Blair, J. T.; Krogh-Jespersen, K. J. Phys. Chem. 1990, 94, 4470. (52) Ando, K.; Kato, S. J. Chem. Phys. 1991, 95, 5966. (53) Cichos, F.; Brown, R.; Bopp, P. A. J. Chem. Phys. 2001, 114, 6834. (54) Arnett, E. M.; McKelvey, D. R. Solvent isotope effect on thermodynamics of nonreacting solutes. In Solute-SolVent Interactions; Coetzee, J. F., Ritchie, C. D., Eds.; Marcel Dekker: New York, 1969; p 343. (55) Lee, S. H.; Rossky, P. J. J. Chem. Phys. 1994, 100, 3334. (56) Beveridge, D. L.; Mezei, M.; Mehrotra, P. K.; Marchese, F. T.; Ravi-Shanker, G.; Vasu, T.; Swaminathan, S. AdV. Chem. Ser. 1983, 204, 297. (57) Ladanyi, B. M.; Skaf, M. S. Annu. ReV. Phys. Chem. 1993, 44, 335. (58) Faeder, J.; Ladanyi, B. M. J. Phys. Chem. B 2000, 104, 1033. (59) Christensen, M.; Pedersen, J. B. J. Chem. Phys. 2003, 119, 5171. (60) Liu, P.; Harder, E.; Berne, B. J. J. Phys. Chem. B 2004, 108, 6595. (61) Svishchev, I. M.; Kusalik, P. G. J. Phys. Chem. 1994, 98, 728. (62) Delbuono, G. S.; Rossky, P. J.; Schnitker, J. J. Chem. Phys. 1991, 95, 3728. (63) Guillot, B.; Guissani, Y. Fluid Phase Equilib. 1998, 151, 19. (64) Guillot, B.; Guissani, Y. J. Chem. Phys. 1998, 108, 10162. (65) Zharkikh, A. A.; Lyashchenko, A. K.; Kharkin, V. S.; Goncharov, V. V.; Lileev, A. S. Zh. Fiz. Khim. 1991, 65, 553.

Zirconia Interface: Molecular

Recommend Documents