MM and Classical Molecular Dynamics Study of [Ru

May 12, 2009 - The solvation of the [Ru(bpy)3]2+ (bpy = 2,2′-bipyridyl) cation was studied by molecular dynamics simulations, using both a hybrid QM...
4 downloads 0 Views 342KB Size
Download PDF
J. Phys. Chem. B 2009, 113, 7737–7744

7737

Combined QM/MM and Classical Molecular Dynamics Study of [Ru(bpy)3]2+ in Water Marc-Etienne Moret,† Ivano Tavernelli, and Ursula Rothlisberger* Laboratoire de chimie et biochimie computationelles, Ecole Polytechnique Fe´de´rale de Lausanne, CH-1015 Lausanne, Switzerland ReceiVed: January 7, 2009; ReVised Manuscript ReceiVed: April 14, 2009

The solvation of the [Ru(bpy)3]2+ (bpy ) 2,2′-bipyridyl) cation was studied by molecular dynamics simulations, using both a hybrid QM/MM Car-Parrinello scheme and a classical force field. The trajectory analysis reveals that the first solvation sphere consists of chains of H-bonded water molecules intercalating between the bpy ligands, which in turns induces an organization of the second solvation sphere. Comparison with a simple charged sphere model shows the influence of the geometry of the cation and its ligands on the surrounding water. A water residence time in the first solvation shell of [Ru(bpy)3]2+ of 12.4 ps is computed, which is significantly larger than the value of 4.5 ps found in bulk water. On the other hand, due to the one-dimensional arrangement of hydrogen bonds, the dipole of the water molecules can easily reorient, providing a way to explain the fast solvent dependent relaxation dynamics observed following photoexcitation. 1. Introduction Polypyridine complexes of ruthenium(II) have been of major interest to the scientific community due to their specific photophysical properties. They exhibit in general relatively lowlying metal-to-ligand charge-transfer (MLCT) states that can readily undergo electron transfer to a suitable acceptor. In particular, they have found applications as photosensitizers for titanium dioxide based solar cells1,2 and as luminescent DNA intercalators.3 The photophysics of ruthenium tris-2,2′-bipyridine, the prototype of this family of complexes, has been extensively studied by spectroscopic4-7 and theoretical8-12 methods. As can be expected from the fact that the MLCT excitation induces an important charge redistribution, the photophysical properties of the [Ru(bpy)3]2+ cation are strongly influenced by the solvent. For example, the rate of localization of the unpaired electron after excitation has been shown to be dependent on the chain length of the nitrile solvent used.13 When water is the solvent, several experimental observations suggest a tight association of water with the complex. The lifetime of the first triplet MLCT (3MLCT) state is much longer in D2O (1.0 µs) than in H2O (0.6 µs). A strong solvent isotope effect on the resonance Raman intensities is observed in water, but not in acetonitrile.14 Finally, the Raman spectrum of [Ru(bpy)3]2+ in zeolites is significantly modified by the presence of water although the size of the cavities barely exceed the size of the cation.15 Thus, an accurate description of the photochemical behavior of [Ru(bpy)3]2+ (bpy ) 2,2′-bipyridyl) requires a thorough understanding of its solvation dynamics. Car-Parrinello molecular dynamics16 is an efficient, density functional theory (DFT) based method for the time-resolved simulation of fairly large systems. Inclusion of a large number of solvent molecules, however, often results in systems that are too large for first-principles simulations. One way to overcome this limitation is to use a hybrid quantum mechanics/molecular * To whom correspondence should be addressed. E-mail: [email protected]. http://lcbcpc21.epfl.ch. † Present address: Laboratorium fu¨r Organische Chemie, ETH Zu¨rich, Wolfgang-Pauli-Str. 10, CH-8093 Zu¨rich, Switzerland.

mechanics (QM/MM) scheme,17,18 where the solute is treated within DFT, while the solvent is described by classical molecular mechanics. For a system of the size considered herein, the time scale accessible for QM/MM simulations is of a few picoseconds. As an accurate statistical analysis of the solvent structure requires longer simulations, we turned to purely classical simulations for that purpose. There are precedents in the literature for molecular mechanics (MM) treatments of ruthenium polypyridyl complexes. Following the development of force-field parameters by the group of Norrby,11 MM methods have been applied to the prediction of structures in the solid state19 and in vacuum,20 as well as to the study of the DNA intercalation of ruthenium polypyridyl compounds.21,22 In this work, the solvation structure and dynamics of ruthenium tris-2,2′-bipyridine in aqueous solution have been investigated by means of DFT and TDDFT based QM/MM simulations in both ground and excited singlet MLCT states. Accurate statistics have been obtained by running purely classical simulations, which were validated by comparison with the QM/MM trajectory. 2. Computational Methods 2.1. DFT and TDDFT Calculations. All ground state (DFT) and excited state (TDDFT) calculations were performed using the Becke exchange functional23 together with the Perdew contribution24 for correlation. A plane wave (PW) basis set with a cutoff of 75 Ry was used. The core electrons were represented by Troullier-Martins25 norm-conserving nonlocal pseudopotentials, in which scalar relativistic effects are included, but no spin-orbit coupling is taken into account. The one electron wave functions were optimized using a direct inversion of the iterative subspace (DIIS)26 algorithm with a wavefunction gradient threshold of 10-5 au; when the latter did not converge, a preconditioned conjugated gradient (PCG) based algorithm was used instead. For the ground state geometry optimization, a DIIS algorithm was used, with a convergence criteria of 5 × 10-4 au for the residual nuclear forces. The ab initio molecular dynamics simulations were performed in the NVT ensemble using the

10.1021/jp900147r CCC: $40.75  2009 American Chemical Society Published on Web 05/12/2009

7738

J. Phys. Chem. B, Vol. 113, No. 22, 2009

Moret et al.

TABLE 1: Parameters Used for Ruthenium in Classical Simulationsa atoms Ru Ru-NC Ru-NC-CA NC-Ru-NC (cis) NC-Ru-NC (trans) Ru-NC-CA-CA Ru-NC-CA-CA Ru-NC-CA-HA NC-Ru-NC-CA (cis)

a

TABLE 2: Charges Used in the MM Simulations of [Ru(bpy)3]2+ in Water

parameters

ref

RvdW ) 2.963 Å ε ) 0.56 kcal · mol-1 req ) 2.081 Å Kr ) 268 kcal · mol-1 · Å-1 θeq ) 123.5° Kθ ) 103.6 kcal · mol-1 · rad-1 θeq ) 91.1° Kθ ) 81.25 kcal · mol-1 · rad-1 θeq ) 180.0° Kθ ) 24.49 kcal · mol-1 · rad-1 (Vn)/(2) ) 1.21 kcal · mol-1 γ ) 0° n)-2 Vn/2 ) 1.41 kcal · mol-1 γ ) 180° n)1 Vn/2 ) 5.27 kcal · mol-1 γ ) 180° n)2 Vn/2 ) 0.25 kcal · mol-1 γ ) 180° n)4

47 11 11 11 11 11

atoma

charge

atoma

charge

Ru N C2 C3 C4 C5

1.465 -0.246 0.041 0.008 0.032 0.018

C6 H3 H4 H5 H6

0.026 0.054 0.067 0.066 0.024

a The atomic numbering for the bipyridine ligand is given in the structure.

11

TABLE 3: Selected Geometrical Features of the Geometry Optimized and the X-ray Structures of [Ru(bpy)3]2+ 11

GO structure

11

The symbols refer to the standard AMBER nomenclature.44

Car-Parrinello scheme16 with a time step of 4 au (∼0.097 fs) and a fictitious electronic mass of 600 au. A constant temperature was imposed by a Nose´-Hoover27,28 thermostat at 300 K with a coupling frequency of 2000 cm-1. TDDFT energies and forces were computed according to refs 29 and 30 within the Tamm-Dancoff approximation (TDA) and the adiabatic local density approximation for the TDDFT kernel (ALDA). A time integration step of 4 au was used for both ground state and TDDFT MD simulations within the Born-Oppenheimer (BO) approximation. All electronic gradients were converged to a maximum deviation of 10-6 au. Oscillator strengths were computed according to Bernasconi et al.31 All DFT and TDDFT calculations were performed with the CPMD code.32 2.2. Classical Simulations. Purely classical (MM) simulations were run using the AMBER833 package. The system contained a [Ru(bpy)3]2+ ion, two Cl- counterions and 3298 TIP3P34 water molecules. As a starting configuration for the ruthenium complex we used the DFT-optimized gas phase geometry. Atomic charges were computed with the Hirshfeld method35 and are listed in Table 2. The complex was represented using parameters from the AMBER/parm99 force field, supplemented by the parameters developed by Norrby et al.11 for ruthenium polypyridyl compounds, to which appropriate unit corrections were applied. The additional parameters used in the simulation are listed in Table 1. The atom types for the organic part were assigned to match as close as possible the chemical properties of each individual atom. According to the AMBER nomenclature, the CA (aromatic carbon) type was assigned to the carbon atoms of the bpy ligand, the NC (sp2 N in 6 membered rings) type to nitrogens and the HA (aromatic hydrogen) to hydrogens. The ability of these parameters to model the compound of interest was tested by comparing the equilibrium distributions of key structural properties with the ones obtained from the QM/MM simulation. For comparison, we also ran a simulation in which the Ru complex was substituted by a charged sphere of approximatively

Ru-N C2-C2′ C2-N C2-C3 C3-C4 C4-C5 C5-C6 C6-N same ligand trans N-C2-C2′-N N-C2-C6-Rua a

X-ray structure48

Bond Lengths (Å) 2.075 2.060 1.467 1.470 1.373 1.364 1.400 1.376 1.391 1.377 1.396 1.361 1.389 1.367 1.355 1.357

% relative deviation 0.7 -0.2 0.7 1.6 1.0 2.6 1.6 -0.1

N-Ru-N Angles (deg) 78.4 78.8 173.1 172.3 Dihedrals (deg) 1.73 -5.56 0.86 2.17

Improper (out-of-plane) angle.

TABLE 4: Mean Values for Selected Geometrical Features with Corresponding Standard Deviations (σ) from the QM/MM Simulation, and Values Obtained for the Optimized Geometry mean value Ru-N C2-C2′ same ligand trans N-C2-C2′-N N-C2-C6-Rua a

σ

Bond Lengths (Å) 2.077 0.046 1.468 0.021 N-Ru-N Angles (deg) 78.4 1.6 172.6 2.8 Dihedrals (deg) 1.5 6.2 1.1 5.1

optimized geom 2.075 1.467 78.4 173.1 1.7 0.86

Improper (out-of-plane) angle.

equivalent radius. The sphere was represented by a new atom type, with a charge of +2, a radius of 5.2 Å, a mass of 569.634 au, and a van der Waals energy well depth of 0.0257 kcal/mol. This pseudoatom was solvated with 3024 TIP3P water molecules, and two chloride ions were added. This system was evolved for 1.1 ns at a constant pressure and temperature of 1 atm and 298 K, respectively. 2.3. QM/MM Setup. We applied a QM/MM molecular dynamics technique based on Car-Parrinello molecular dynamics16 that has been demonstrated to provide an accurate description of solvation properties of neutral as well as charged systems.17,18,36-38 The steric interaction between the QM and

Molecular Dynamics Study of [Ru(bpy)3]2+ in Water

J. Phys. Chem. B, Vol. 113, No. 22, 2009 7739

the MM part is modeled by a Lennard-Jones potential as described by the AMBER force-field, while the electrostatic interactions with the neighboring MM atoms are taken into account by the energy term

Eel )

∑ qi ∫ dr F(r) Vi(|r - ri|)

(1)

i

where qi is the charge of MM atom i, ri is its position, F ) Fel + Fion is the total (ionic and electronic) charge density of the quantum system, and Vi is the Coulomb potential of atom i modified at short-range to avoid spurious overpolarization effects.17 The functional form of Vi is given by

Vi(r) )

n rc,i - rn n+1 rc,i - rn+1

(2)

where n is 4 and rc,i is equal to 0.32 Å for hydrogen and to 0.72 Å for oxygen. The electrostatic interaction between the QM system and distant MM atoms is modeled by coupling of the multipole moments of the quantum charge distribution with the classical point charges.17 The classical electrostatic potential polarizes the charge density of the QM system and thus affects its electronic properties. The QM subsystem (Ru complex) was treated at the DFT/BP level as described above, while water (TIP3P) and chloride ions were treated classically using the AMBER force field. The inherent periodicity in the plane wave calculations was circumvented by solving Poisson’s equation for nonperiodic boundary conditions,39 while periodic boundary conditions were retained for the classical solvent box. The electrostatic interaction is taken into account by the P3M method.40 A starting structure for the QM/MM simulation was obtained from a classical 1.1 ns simulation with a fixed geometry for the ruthenium complex performed at a constant pressure of 1 atm and a temperature of 298 K with 3024 TIP3P water molecules. In this starting frame, the two chloride anions were at distances of 23 and 24 Å from the ruthenium ion, and at 18 and 19 Å from the closest atom of the complex, and no ion pairing was observed throughout the QM/MM simulation. The system was first thermalized at increasingly higher temperatures from 100 to 300 K in three intervals of 100 K with fixed solute geometry. Subsequently, all constraints were removed and the system was equilibrated for 2.5 ps at 300 K. All these preparatory relaxation and heating steps were performed using Born-Oppenheimer dynamics, while production dynamics in the ground state were produced using the Car-Parrinello electronic propagation scheme. Born-Oppenheimer dynamics was also used in the TDDFT-QM/MM calculations18,41 with a time step of 5 au. All QM/MM simulation were performed in the NPT ensemble at 1 atm and 300 K. 3. Results and Discussion 3.1. Optimized Geometry in Vacuo. As shown in Table 3, the geometry optimized complex shows few structural differences from the starting crystal structure, in agreement with the observation by Buchs and Daul.9 Concerning the intraligand NCCN dihedral angle, our geometry optimized structure shows an almost flat bipyridine ligand, while the crystal structure shows a N-C-C-N dihedral angle of ∼-5.6°. Furthermore, this dihedral angle has changed sign during the geometry optimization, which shows that the small but nonzero angle of ∼1.7° is not due to a too loose convergence criterion.

The observed small deviation from the crystal structure may be assigned to the loss of crystal packing effects in our calculation, where the molecule is isolated. Buchs and Daul calculated a dihedral angle of ∼7.1°, which might be thought to infirm this hypothesis. The discrepancy between our results and the previous ones reported in ref 9 may be due to a bias introduced by the choice of the starting geometry for the isolated ligand (cis configuration with a clear nonplanarity42) used to model the complex, and by the accuracy chosen for the convergence of the atomic gradients (1 × 10-2 in ref 9 vs 5 × 10-4 hartree/Å in this work). In addition, as observed in our QM/MM simulation, the NCCN dihedral angle shows large amplitude fluctuations at room temperature; we measured an average value of 1.5° with a standard deviation of about 6°. 3.2. Classical vs QM/MM Structural and Dynamical Properties. The hybrid Car-Parrinello/MM simulation of the complex in its electronic ground state dissolved in water was run at 300 K for 2.9 ps. The mean values and standard deviations of some selected geometrical features are given in Table 4. The mean values computed for all parameters considered in this study are close to the ones measured for the optimized gas phase geometry, which means that solvation in water does not induce large structural changes of the solute. The coordination geometry around the ruthenium ion is found to be quite rigid: the Ru-N bond length has a standard deviation of 0.046 Å (∼2.2%), and the trans N-Ru-N angle has one of 2.8°. Another interesting fact is that the mean value of the characteristic torsional angle of the ligand (the N-C2-C2′-N dihedral angle) is not zero, but very close to that found in the optimized geometry (1.7°). The flexibility of this dihedral angle might seem surprising if we consider both the loss in π conjugation and the nonideal orientation of the nitrogen lone pair induced by a torsion. However, those two effects are partially compensated by the steric repulsion between the 3 and 3′ hydrogen atoms of the bpy ligand (2.09 Å distance in the optimized geometry), which causes the free ligand not to be planar in its cis conformation.42 To extend the simulated time scale and to get a clearer picture of the long-time solvation dynamics of [Ru(bpy)3]2+ in aqueous solution, a ∼1 ns long classical molecular dynamics simulation was performed using the AMBER33 package. The geometry around the ruthenium center was modeled using the bonding parameters (equilibrium bond lengths and angles as well as the corresponding force constants and dihedral potentials) developed by Norrby et al.11 for the force field MM3*, used in AMBER with appropriate unit corrections. The functional form of MM3*43 is similar to that of AMBER,44 except for the inclusion of higher order terms for bond and angle stretching in the former. Because this may in principle lead to some transferability issues, the used parametrization was validated by comparison with the QM/MM trajectory. To check the quality of our parametrization of the solute, the equilibrium distribution of some key structural parameters were compared with the ones obtained from the QM/MM simulation. This comparison is shown in Figure 1. Most characteristic geometrical features are indeed correctly represented by our classical model. The main exception concerns the bipyridine bridging carbon-carbon bond, which has clearly a shorter equilibrium length and a smaller force constant in the classical dynamics representation. This can easily be explained by the fact that in our setup the same atom type was assigned to every carbon atom, which means that this specific bond is treated as if it would belong to an aromatic cycle. This problem could have been circumvented by assigning new parameters to

7740

J. Phys. Chem. B, Vol. 113, No. 22, 2009

Figure 1. Comparison of the thermal distributions for selected distances (Å), angles (deg), and dihedral angles (deg) obtained from the MM (black) and the QM/MM (red) simulations.

the bridging carbon atoms, but we expect that a change of 0.05 Å in this particular bond length would not strongly affect the solvation dynamics. We also observed a difference in the distribution of the N-Ru-N angle in trans position. This can be understood considering that, according to the force field parameter set used to perform these calculations, all cis N-Ru-N angles are assigned a single force constant. This may induce an error in the geometrical arrangement of the six nitrogens around the metal. As can be noticed from Figure 1, this effect is imperceptible in the distribution of angles involving two nitrogen atoms belonging to the same ligand, because this angle is constrained by the intrinsic ligand geometry. Furthermore, we assume that the influence of a small deviation in the coordination geometry resulting in a difference of about 5° in this angle does not have a strong influence on the solvent structure around the complex. 3.3. Analysis of the Solvation Structure. The radial distribution functions g(r) of water oxygen and hydrogen atoms around the ruthenium have been computed together with the corresponding distance dependent coordination number N(r) for both the QM/MM and MM trajectories (Figure 2). The main difference between the g(r) functions computed from the two trajectories is the stronger noise observed for the QM/MM data, which is due to the shorter accessible time scale. Since the QM/MM trajectory duration (2.9 ps) is much smaller than the water residence time in the first solvation shell (12.4 ps, vide infra), the sampling of the phase space is insufficient. Thus, while the polarizability introduced by the quantum treatment of the solute makes the description of the system in principle more accurate, the shorter accessible time scale renders the Car-Parrinello QM/MM method less appropriate for an accurate description of the time-averaged solvent structure. However, we observe that the overall structure of the radial distribution functions is the same for both methods. In the following paragraphs, the solvent structure around the [Ru(bpy)3]2+ cation will be discussed on the basis of the MM trajectory analysis.

Moret et al.

Figure 2. Radial distribution functions g(r) and running coordination numbers N(r) (dashed lines) of the water oxygen (black) and hydrogen (red) atoms with respect to the ruthenium center, computed from the QM/MM (A) and MM (B) trajectories. (C) g(r) and N(r) computed for a pseudo atom of charge + 2 and van der Waals radius 5.2 Å.

Figure 3. Plot of the ratio NH(r)/NO(r) as a function of the distance from the ruthenium ion, computed from the QM/MM (dashed line) and MM (solid line) trajectories, as well as for the sphere model (dotted line).

The Ru-O distribution exhibits a relatively broad peak around 5.6 Å, followed by a minimum at about 6.5 Å that corresponds to a first solvent shell made of ∼15 water molecules. These molecules are closer to the ruthenium center than the most external hydrogens of the bipyridine moieties and thus can be interpreted as intercalation between the ligands. The first peak of the Ru-H g(r) is slightly shifted further away from the center, and also a little broader. This reflects a preferential orientation of the water dipole toward the metal. A more precise picture of the orientation of the water molecules is obtained by plotting the ratio nHO(r) ) NH(r)/NO(r) (where NH(r) and NO(r) are the running coordination numbers of hydrogen and oxygen atoms, respectively) as a function of the distance to the ruthenium center, as shown in Figure 3. This quantity represents the ratio between the number of hydrogen and oxygen atoms contained in a sphere of radius r centered at the position of the ruthenium atom. Because of the smaller van der Waals radius of the hydrogen, nHO(r) has a high value (>2) at very short r, even though both values NH(r) and NO(r) are smaller than 1 in this range. As r increases, nHO(r) drops below 2, reflecting the orientation of the oxygen atoms toward the positive charge of the ruthenium complex, and then tends

Molecular Dynamics Study of [Ru(bpy)3]2+ in Water

Figure 4. Snapshot from the MM simulation showing the structure of the first water shell. The atoms are represented by spheres of standard van der Waals radii, with the following color code: white, H; red, O; cyan, C; blue, N; green, Ru. The image represents all atoms in a radius of 6.9 Å from the ruthenium with the exception a few hydrogen atoms from the second shell that have been removed for clarity.

asymptotically to the bulk value of 2. However, in the case of [Ru(bpy)3]2+, nHO(r) always stays above 1, which means that within a given distance from the ruthenium ion there are always at least as many hydrogen as oxygen atoms. This fact can be explained by the formation of a linear chain of water molecules in the first solvation shell linked together by hydrogen bonds and with one hydrogen per molecule pointing outside. In fact, fragments of such a tennis ball like arrangement of water molecules in the groove between the ligands are observed on virtually all frames obtained from our simulations (an example is given in Figure 4). The geometrical distribution of the first solvation shell around the complex is represented in Figure 5. This plot is obtained by sampling coordinates of the oxygen atoms within a radius of 6.5 Å of the ruthenium center along the trajectory according to an internal orthonormal coordinates system defined as shown in Figure 6. The xy plane is defined by the centers of the three C2-C2′ bonds of the bpy ligands, with the origin at their barycenter and the unitary vector b x pointing away from one of them. In the frozen geometry, the origin corresponds exactly to the Ru atom, and the b z axis to the C3 axis of the complex. This plot confirms that the first peak in the radial distribution function corresponds to water molecules intercalating between two or three ligands. A more detailed picture of the structure of the first solvation shell was obtained by considering the dependency of the connectivity of the water molecules on the angle θ, which is defined as the angle between the Ru-O vector and the z axis. Comparison of the computed distribution of the water molecules (Figure 7, bars) with that corresponding to a spherical symmetry (Figure 7, dashed line) reveals three highly populated regions. Two of them (θ < 30° and θ > 150°) are caused by the cavities formed by the three bpy units around the D3 axis, while the third one (66° < θ < 114°) corresponds to water molecules intercalating betwen two ligands. For each value of θ, we computed the average number of hydrogen bonds45 a water molecule donates to (circles) and accepts from (squares) another water molecule from the first (red) or second (blue) solvation

J. Phys. Chem. B, Vol. 113, No. 22, 2009 7741 shell (Figure 7, top). We find that a higher fraction of the hydrogen bonds are connecting two water molecules from the first shell around θ ) 90°, which we interpret as a consequence of the formation of one-dimensional chains in the interligand groove. This interpretation is supported by the fact that the average number of hydrogen bonds per water molecule from the first to the second shell is high (ca. 0.8) compared to that of hydrogen bonds from the second to the first shell (ca. 0.45). As expected from the low polarity of the C-H bonds, virtually no hydrogen bonds were found from the bipyridine ligand to water molecules. The peculiar structure of the first shell is expected to induce some organization in the second shell. In fact, the g(r) functions plotted in Figure 2B exhibit a quite well-defined second shell peak around 8 Å for both hydrogen and oxygen, and even a smooth but still clearly identifiable third peak. The water molecules in these shells show very little tendency for orientation with respect to the ruthenium center, which means that these structures are likely not due to long-range electrostatic effects only. In one of the simplest solvation models one could think of, [Ru(bpy)3]2+ can be represented by a charged sphere with a radius comparable to the one of the solute cavity. To test the validity of such a model and to distinguish between electrostatic and geometric effects on the solvation structure, a further simulation was performed under the same conditions but with the complex replaced by a pseudoatom with a van der Waals radius of 5.2 Å and an energy well depth of 0.0257 kcal/mol. The chosen radius corresponds to the Ru-H5 distance, which was thought to reasonably model the size of the ion of interest. This choice is otherwise arbitrary, but as we are interested mostly in qualitative observations of the solvent structure, slightly different choices of the radius will not affect our results. The same analysis as for the complex was performed, and the results are presented in Figure 2C. The solvation picture for such a model is quite different from the one we get from the explicit solute simulation. The most striking differences are related to the magnitude of the first shell peak and to the strong orientation of the water molecules, which appears as an outward shift of ∼0.5 Å of the hydrogen peak with respect to the oxygen peak. In the case of this model system, the first solvation shell thus contains about 35 well oriented water molecules. The stronger orientation of the water molecules in the first shell is evident from the plot in Figure 3; the nHO(r) ratio exhibits a minimum as low as 0.52 in the first shell region. The second shell is poorly definedsthe corresponding peaks in the g(r) functions are broadsand the solvent shows almost no structure beyond 10 Å from the center of the pseudo-atom. From these observations one can conclude that the specific solvent structure around [Ru(bpy)3]2+ in aqueous solution is not only due to electrostatic interactions with the +2 charge of the ion but also due to the water intercalation and the contact with the apolar surfaces of the ligands. Considering the fact that the first solvation shell consists mainly of intercalating water molecules (Figure 5), one may suggest that the sphere model would better represent the complex ion and the first solvation shell. However, a qualitative comparison between the second peak of Figure 2B and the first one of Figure 2C shows that this model is of limited accuracy: the charged-sphere model induces a strong orientation of the water molecules while the molecules in the second shell of [Ru(bpy)3]2+ are almost not oriented. 3.4. TDDFT Calculations. The peculiar geometrical arrangement of the first solvation shell of the complex together with its dynamical properties may play an important role in the

7742

J. Phys. Chem. B, Vol. 113, No. 22, 2009

Moret et al.

Figure 5. Position of the oxygen atoms in a radius of 6.5 Å from ruthenium sampled along a 1 ns classical trajectory. The positions are given with respect to the internally defined orthonormal coordinates system described in Figure 6. Left: projected along the y axis. Right: projection along the z axis. A sphere with a radius of 4.5 Å is represented for a better visualization.

Figure 6. Internal coordinate system definition.

Figure 8. Relaxation dynamics of the [Ru(bpy)3]2+ complex in water upon excitation into the first optically active singlet MLCT. The four snapshots are taken from a TDDFT/QM/MM trajectory at 5 (upper left), 20 (upper right), 40 (lower left), and 60 (lower right) fs after excitation. In green are shown the two water molecules centered on top of the bipyridine rings, which undergo a major and ultrafast reorientation following photoexcitation.

Figure 7. Structure of the first solvation shell as a function of the angle θ, defined as the angle between the Ru-O vector and the z axis defined in Figure 6. Bottom: average number of water molecules per 6° conical slice. The dashed line represents an ideal spherical distribution. Top: average number of hydrogen bonds donated to (circles) or accepted by (squares) each water molecule in the first solvation shell. Key: red, H-bonds within the first shell; blue, H-bonds between the first and the second shell.

relaxation processes following photoexcitation. Upon excitation, an electron is transferred from the metal to the ligands leading to a singlet MLCT complex. This induces an important change

in the electric field felt by the surrounding solvent molecules, which will reorient to minimize the interaction with the complex. To investigate this process, we have performed three independent TDDFT excited state MD calculations using the same QM/ MM setup employed in the ground state calculations. The system was first vertically excited into the singlet MLCT state with the largest oscillation strength and then propagated for 250 fs using the TDDFT BO MD scheme implemented in CPMD.18 Figure 8 shows four snapshots taken at times 5, 20, 40, and 60 fs after excitation along one of the simulated excited state trajectories. For clarity, only the first and part of the second solvation shells are shown. The complex is oriented with the ligand carrying the largest fraction of the extra charge pointing toward the reader. For two water molecules in the first solvation shell of the charged ligand we observe a very fast reorientation consisting in a rotation of the “free” OH bond toward the center of the underneath pyridine ring. This ultrafast relaxation process is only possible due to the quasi one-dimensional arrangement of the water molecules in the grooves between the ligands.

Molecular Dynamics Study of [Ru(bpy)3]2+ in Water

J. Phys. Chem. B, Vol. 113, No. 22, 2009 7743

TABLE 5: Residence Time, τs (ps), for the Water Molecules in the First Solvation Shell of [Ru(bpy)3]2+, the Charged Sphere Model, and in Bulk Water τs

center [Ru(bpy)3] sphere2+ H 2O

2+

References and Notes

12.4 8.4 4.5

3.5. Water Residence Times. The dynamic properties of the first coordination shell were investigated by the method described by Impey et al.46 In this scheme, a function nion is defined as Nt

t* nion (t)

∑ ∑ Pjt (tn,t)

1 ) Nt n)1

*

Acknowledgment. We acknowledge useful discussions with Majed Chergui and Tony Vlcek.

(3)

j

where Pt* j (tn,t) is equal to 1 if molecule j lies within the first solvation shell of the ion at both tn and tn + t without leaving it for more than t*, and to 0 in all other cases. nion(0) is by definition equal to the coordination number of the solute ion. Except at very short times, the function nion(t) can be accurately fitted by a exponential decay function of type nion(t) = nion(0) s s )), of which the time constant τion is the residence exp(-(t)/(τion time in the considered shell. The parameter t* accounts for molecules which temporarily leave the shell without properly entering the bulk and come back shortly after. We applied this analysis to the first solvation shell of [Ru(bpy)3]2+, defined as water molecules of which the oxygen atom is located within 6.3 Å of the ruthenium center. The value s 2+ naturally increases with increasing values of the t* of τ[Ru] parameter. The increase was found to become less important for t* g 2.25 ps, and this value was chosen for subsequent analysis. (Impey et al.46 use a similar value of t* ) 2 ps.) The obtained residence times are presented in Table 5. As can be seen from these values, the water chain-like structure around [Ru(bpy)3]2+ is short-lived, but water is still retained about 3 times longer in the first solvation sphere than in bulk water. Interestingly, the charged sphere model exhibits faster water exchange than the complex, showing that electrostatic effects alone cannot account for the observed exchange rate. 4. Conclusions The aqueous solvation of the [Ru(bpy)3]2+ cation has been investigated by means of QM/MM and purely classical molecular dynamics simulations. The existence of a first hydration shell consisting of ca. 15 water molecules intercalated between the bpy ligands and a structured second shell was demonstrated. The first shell was observed to form chains of hydrogen bonded molecules placed along the grooves in between the ligands. However, this structure is highly dynamic, as the characteristic water-exchange time could be evaluated to be ∼12 ps. These observations shed some new light onto the “tight association” of the water molecules with the [Ru(bpy)3]2+ cation that is evident from experiment. The undercoordination of these water molecules is responsible for their interesting dynamical properties, which may play an important role in determining the photophysical properties of these type of complexes. Among others, this one-dimensional arrangement of hydrogen bonds allows for a fast reorientation of the water dipoles, providing a way for an efficient and fast stabilization of the MLCT state induced by photoexcitation. These results will hopefully contribute to a better understanding of the fascinating photochemistry of ruthenium polypyridine complexes.

(1) Gra¨tzel, M. Nature 2001, 414, 338–344. (2) Hagfeldt, A.; Gra¨tzel, M. Acc. Chem. Res. 2000, 33 (5), 269–277. (3) Metcalfe, C.; Thomas, J. A. Chem. Soc. ReV. 2003, 32, 215–224. (4) Saes, M.; Bressler, C.; Abela, R.; Grolimund, D.; Johnson, S. L.; Heimann, P. A.; Chergui, M. Phys. ReV. Lett. 2003, 90 (4), 047403. (5) Damrauer, N. H.; Cerullo, G.; Yeh, A.; Boussie, T. R.; Shank, C. V.; McCusker, J. K. Science 1997, 275, 54–57. (6) Juris, A.; Balzani, V.; Barigelletti, F.; Campagna, S.; Belser, P.; Von Zelewsky, A. Coord. Chem. ReV. 1988, 85–277. (7) Yersin, H.; Humbs, W.; Strasser, J. Coord. Chem. ReV. 1997, 159, 325–358. (8) Daul, C.; Baerends, E. J.; Vernooijs, P. Inorg. Chem. 1994, 33, 3538–3543. (9) Buchs, M.; Daul, C. Chimia 1998, 52, 163–166. (10) Zheng, K. C.; Wang, J.; Peng, W. L.; Liu, X. W.; Yun, F. C. J. Mol. Struct. (THEOCHEM) 2002, 582, 1–9. (11) Brandt, P.; Norrby, T.; Åkermark, B.; Norrby, P.-O. Inorg. Chem. 1998, 37, 4120–4127. (12) Damrauer, N. H.; Weldon, B. T.; McCusker, J. K. J. Phys. Chem. A 1998, 102, 3382–3397. (13) Yeh, A. T.; Shank, C. V.; McCusker, J. K. Science 2000, 289, 935– 938. (14) Adam Webb, M.; Knorr, F. J.; McHale, J. L. J. Raman Spectrosc. 2001 2001, 32, 481–485. (15) Incavo, J. A.; Dutta, P. K. J. Phys. Chem. 1990, 94 (7), 3075– 3081. (16) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55 (22), 2471–2474. (17) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Chem. Phys. 2002, 116 (16), 6941–6947. (18) Moret, M.-E.; Tapavicza, E.; Guidoni, L.; Rohrig, U. F.; Sulpizi, M.; Tavernelli, I.; Rothlisberger, U. Chimia 2005, 59, 493–498. (19) Breu, J.; Domel, H.; Norrby, P.-O. Eur. J. Inorg. Chem. 2000, 2409– 2419. (20) Abrahamsson, M.; Lundqvist, M. J.; Wolpher, H.; Johansson, O.; Eriksson, L.; Bergquist, J.; Rasmussen, T.; Becker, H.-C.; Hammarstro¨m, L.; Norrby, P.-O.; Åkermark, B.; Persson, P. Inorg. Chem. 2008, 47 (9), 3540–3548. (21) Han, D.; Wang, H.; Ren, N. J. Mol. Model. 2004, 10, 216–222. (22) Han, D.; Wang, H.; Ren, N. J. Mol. Struct. 2004, 711, 185–192. (23) Becke, A. D. Phys. ReV. A 1988, 38 (6), 3098–3100. (24) Perdew, J. P. Phys. ReV. B 1986, 33 (12), 8822–8824. (25) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43 (3), 1993–2006. (26) Hutter, J.; Lu¨thi, H. P.; Parrinello, M. Comput. Mater. Sci. 1993, 2, 244–248. (27) Nose´, S. Mol. Phys. 1984, 52 (2), 255–268. (28) Hoover, W. G. Phys. ReV. A 1985, 31 (3), 1695–1697. (29) Hutter, J. J. Chem. Phys. 2003, 118 (9), 3928–3934. (30) Tavernelli, I.; Ro¨hrig, U. F.; Rothlisberger, U. Mol. Phys. 2005, 103, 963–981. (31) Bernasconi, L.; Sprik, M.; Hutter, J. Chem. Phys. Lett. 2004, 394, 141–146. (32) CPMD, IBM Corp. 1990-2004, MPI fu¨r Festko¨rperforschung Stuttgart 1997-2001. (33) Case, D.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Wang, B.; Pearlman, D. A.; Crowley, M.; Brozell, S.; Tsui, V.; Gohlke, H.; Mongan, J.; Hornak, V.; Cui, G.; Beroza, P.; Schafmeister, C.; Caldwell, J. W.; Ross, W. S.; Kollman, P. A. AMBER; 2004. (34) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79 (2), 926–935. (35) Hirshfeld, F. L. Theor. Chim. Acta 1977, 44, 129–138. (36) Gossens, C.; Tavernelli, I.; Rothlisberger, U. Chimia 2005, 59, 81– 84. (37) Gossens, C.; Tavernelli, I.; Rothlisberger, U. J. Chem. Theory Comput. 2007, 3, 1212–1222. (38) Masson, F.; Laino, T.; Tavernelli, I.; Rothlisberger, U.; Hutter, J. J. Am. Chem. Soc. 2008, 130, 3443–3450. (39) Martyna, G. J.; Tuckerman, M. E. J. Chem. Phys. 1999, 110, 2810– 2821. (40) Luty, B. A.; van Gunsteren, W. F. J. Chem. Phys. 1996, 100, 2581– 2587. (41) Cascella, M.; Cuendet, M. A.; Tavernelli, I.; Rothlisberger, U. J. Phys. Chem. B 2007, 111, 10248–10252. (42) Howard, S. T. J. Am. Chem. Soc. 1996, 118, 10269–10274. (43) Allinger, N.; Yuh, Y. H.; Lii, J.-H. J. Am. Chem. Soc. 1989, 111 (23), 8551–8566.

7744

J. Phys. Chem. B, Vol. 113, No. 22, 2009

(44) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz Jr, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117 (19), 5179–5197. (45) A hydrogen bond was defined as an oxygen-hydrogen contact with an O-H distance shorter than 2.30 Å and an O-H-O angle between 150° and 210°. (46) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071–5083.

Moret et al. (47) Adlhart, C.; Chen, P. Angew. Chem., Int. Ed. 2002, 41 (23), 4484– 4487. (48) Rillema, D. P.; Jones, D. S.; Levy, H. A. Chem. Commun. 1979, 849.

JP900147R