J . Phys. Chem. 1988, 92, 1008-1010 One last remark about metal catalysis should be mentioned. The above results suggest that if one were to look for a metal surface to catalyze the MCN but not MBr formation, the structure of the active site on this surface has to be made of as few metal atoms as possible.
Acknowledgment. We thank the Office of Naval Research for financial support. Assistance of Drs. E. L. Chronister and R. J. St. Pierre is greatly acknowledged. M. A. El-Sayed thanks Prof. R. B. Bernstein for useful discussions. A. Eychmuller also thanks the Deutsche Forschungsgemeinschaft for a research scholarship.
Monte Carlo Study of the Dlelectric Saturation in Molecular Solutions. Physical Basis of the C Mode in Electron-Transfer Reactlons Yasuyo Hatano, Faculty of Liberal Arts, Chukyo University, Shouwa- ku, Nagoya 466, Japan
Minoru Saito, Toshiaki Kakitani,* Department of Physics, Nagoya University, Chikusa- ku, Nagoya 464, Japan
and Noboru Mataga* Department of Chemistry, Faculty of Engineering Science, Osaka University, Toyonaka Osaka 560, Japan (Received: October 8, 1987; In Final Form: December 18, 1987)
By a Monte Carlo simulation of polar solution using a spherical hard-core model for molecules, we have demonstrated that a thin dielectric saturation shell exists just around a charged solute molecule. The curvature of the free energy curve as a function of solvent polarization corresponding to this dielectric saturation shell is found to be greatly increased with increase of the charge of the solute molecule. The reorientational solvent motion with this property of dielectric saturation is deemed to be the coordinated solvent mode (c mode) which has played a crucial role in the energy gap law of the electron-transfer reaction.
Introduction Electron transfer (ET) is the most fundamental process in chemistry and biology as well as in physics. Many theoretical and experimental studies have been made so far for the various phases of materials.Ib3 Among them, the ET in polar solution has been a field of current lively investigation.eg The regulation mechanism of the E T rate by means of the free energy difference between the initial and final states, called the energy gap law, was one of the central problems.e15 That is, the existence or nonexistence of the "inverted region" at the large energy gap has been extensively investigated.'*15 In this respect, we have recently proposed a theoretical treatment with a new aspect of the solvent motion that the potential curvature for the reorientation of solvent molecules around a charged solute molecule (reactant or product) will be much larger than that around a neutral solute molecule (1) Devault, D. Q.Reu. Biophysics. 1985, 13, 387. (2) Tunneling in Biological Systems; Chance, B., et al., Eds.; Academic: New York, 1979. (3) Beitz, J. V.; Miller, J. R. J . Chem. Phys. 1979, 71, 4579. (4) Marcus, R. A. J. Chem. Phys. 1956, 24, 966. (5) Marcus, R. A. Annu. Rev. Phys. Chem. 1964, IS, 155. (6) Dogonadze, R. R. Dokl. Chem. 1959, 124, 9. (7) Levich, V. G. Ado. Electrochem. Electrochem. Eng. 1966, 4, 249. (8) Kestner, N. R.; Logan, J.; Jortner, J. J. Phys. Chem. 1974, 78, 21. (9) Ulstrap, J.; Jortner, J. J . Chem. Phys. 1975, 63, 4835. (IO) (a) Rehm, D.; Weller, A. Isr. J . Chem. 1970.8, 259. (b) Ballardini, R.; Varani, G.; Indelli, M. T.; Scandola, F.; Balzani, V. J . Am. Chem. SOC. 1978, 100, 7219. (c) Vogelman, E.; Schreiner, S.;Rauscher, W.; Kramer, H. E. A. Z . Phys. Chem. (Frankfurt am Main) 1976, 101, 321. (1 1) Wasielewski, M. R.; Niemczyk, M. P.; Svec, W. A.; Pewitt, E. B. J . Am. Chem. Soc. 1985, 107, 1080. (12) Miller, J. R.; Calcaterra, L. T.; Closs, Q. L. J . Am. Chem. Soc. 1984, 106, 3047. (13) (a) Mataga, N.; Kanda, Y.; Okada, T. J. Phys. Chem. 1986,90,3880. (b) Mataga, N.;Shioyama, H.; Kanda, Y. J. Phys. Chem. 1987,91, 314. (c) Mataga, N. Acta Phys. Polon. 1987, A l l , 767. (14) (a) Ohno, T.; Yoshjmura, A.; Mataga, N. J. Phys. Chem. 1986,90, 3295. (b) Ohno, T.; Yoshimura, A.; Shioyama, H.; Mataga, N. J. Phys. Chem. 1987, 91, 4365. (15) Kakitani, T.; Mataga, N. J. Phys. Chem. 1986, 90, 993
0022-3654/88/2092-1008$01.50/0
(reactant or product).I6 This solvent mode is called the solvent-coordinated mode (c mode, for brevity). Taking into account the c-mode, we have derived three different energy gap laws with respect to the three types of the ET reaction:'s19 (i) no appreciable inverted region for the charge separation reaction (A-B A p B s + ) , (ii) moderate inverted region for the charge shift reA-B;), and (iii) quite remarkable inverted action (A>-B A-B). region for the charge recombination reaction ( A p B ; All those theoretical predictions were consistent with most of the available experimental data."15 Since the above theoretical treatments were based on the assumption of a large frequency change of the solvent motion around a molecular reactant, it is necessary to prove the validity of its assumption more directly. At this point, it should be noticed that so far as the solvent polarization is linear to an external field, the force constant of the solvent motion around the neutral molecule is the same as that around the charged molecule. This fact indicates that a nonlinear response of the solvent polarization, probably due to the solvent dielectric saturation, must be responsible to a large change of the force constant of the solvent motion. In the following, we conduct the Monte Carlo simulation of polar solution and investigate to what extent the dielectric saturation is realized and the force constant is changed, depending upon the magnitude of a charge of the solute molecule.
-
-
-
Method of Monte Carlo Simulation We adopt the following simple model for a polar solution: A solute molecule has a charge re and its shape is a sphere with a hard-core radius roe The solvent molecule has a point dipole moment p and it has a spherical hard core with a radius a. The molecular ion is always fixed at the center of a large spherical (16) Kakitani, T.; Mataga, N. Chem. Phys. 1985, 93, 381. (17) Kakitani, T.; Mataga, N. J. Phys. Chem. 1985, 89, 8 . (18) Kakitilni, T.; Mataga, N. J . Phys. Chem. 1985,89, 4752. (19) Kakitani, T.; Mataga, N. J . Phys. Chem. 1987, 91, 6277.
0 1988 American Chemical Society
The Journal of Physical Chemistry, Vol. 92, No. 5, 1988 1009
Letters
-DI
a
-1.7-
z= 0
: 1
-3.4 1 L '
0 , Fa
-1.7
z=O.5
Q
'0
0 2
0 4
0 6
0 8
1 0
1 2
1 4
1 6
1 8
lfW
2 0
z
Figure 1. Calculated results of P J p as a function of z for some shells. The number in the figure denotes the shell number. The shell width is 1.1 A.
vessel with a radius R. Nosolvent molecules are put into this vessel and are allowed to move under a restriction that the centers of solvent molecules are inside the spherical vessel of the radius R. Before going into the actual Monte Carlo simulation, we first mix those solvent molecules considerably (ca. 1000 steps), neglecting the dipole-dipole and dipole-charge interactions. Then, the position and orientation of n solvent molecules inside the boundary shell between the radii of R and R - 3.5 A are fixed. For the residual N (=No - n) solvent molecules, Monte Carlo simulation is conducted by restoring all the interactions. The n boundary solvent molecules are introduced as a check of the validity of our simulation for a finite size of solution. It should be noticed that the periodic boundary condition with use of the Ewald-Kornfeld method20 is not strictly applicable becai se the system has a net charge and that the actual calculation by this method has been done so far by adopting a certain truncation of the long-range interaction.21 We adopt the following typical parameter values in molecular solutions: a = r, = 2.2 A, R = 25 %, L./
= 2 D, N o = 6 4 4
N = 360,
9
= 0.45
T = 300 K, MCS = 104
(1)
where 7 and MCS denote the space-filling ratio of the solvent and the total number of Monte Carlo steps, respectively.
Results of Monte Carlo Simulation Our simple model system seems to go to equilibrium very quickly when z 2 0.5: The total energy goes into the level of thermal fluctuation in the first 500 M C steps. In the present analysis, we neglect the initial 1000 Monte Carlo steps and use remaining 9000 steps for z 1 0.5, but we neglect the initial 7000 steps and use remaining 3000 steps for z < 0.5. W e divide a space inside the spherical vessel into some shells with width 1, starting from a radius ro + a and call the respective shell from the inside 1st shell, 2nd shell, and so on for convenience. We calculate the polarization in the radial direction normalized to one solvent molecule for each shell and take an ensemble average. The polarization calculated in this way is represented by PI. We plot P r / p as a function of z in Figure 1. The shell width 1 was chosen to be 1.1 A. It should be noticed that the magnitude of z stands for the strength of the external field to the solvent system. It is clearly seen that the dielectric saturation begins to occur a t z = 0.5 and a considerable saturation takes place at z = 1.0 and the saturation is nearly complete at z = 2.0 in the first shell. In contrast to this, no appreciable dielectric saturation can be seen (20) (a) Ewald, P. P. Ann. Phys. 1921, 64, 253. (b) Kornfeld, H. Z . Phys. 1924, 22, 21. (21) Engstrom, S.; Jonsson, B.; Impey, R. W. J . Chem. Phys. 1984,80, 5481.
" -1.7
z=l
.o
-3.4
-3.4 -1.0
'
-0.8'
-6.6 '-0.4 ' - 0 2
'
'0 P,/
'
62
'
0:4
0:s
'
0s:
' 1.0
CI
Figure 2. Calculated free energy curves as a function of P r / p for the first shell with variation of z . The unit of the free energy is kcal/mol.
for the other shells at least until z = 1. An unexpected result is obtained for the 2nd shell that the (Pr/p)-z curve for the second shell is concave at z = 0.6-1.0. A similar behavior is also weakly seen for the 3rd shell. The behavior of the 2nd and 3rd shells may be qualitatively understood as follows: As the value of z increases more than 0.5, the dielectric saturation in the first shell proceeds promptly. Then, the electric field passing through the first shell becomes much stronger than that in the case of the nonsaturation, because the screening effect is weakened. This stronger field causes the greater polarization in the 2nd and 3rd shells where the linear-response still holds. The concave behavior of the curve for the 2nd and 3rd shells ceases at the larger value of z , which must reflect the beginning of a weak dielectric saturation in these shells. We can confirm this fact later in the calculation of curvature. Next, we calculate the free energy curve as a function of the solvent polarization. For this purpose, we make a histogramf for the population of Pr/p in each shell, which resides in the region Pr/K to (Pr/p)-!- ( A P J K ) , in the total Monte Carlo steps ( A P r / p being chosen as 11400). Since ff(Pr/p) is proportional to the probability distribution that Pr/p is realized, it should be expressed by the Boltzmann distribution Y(PI/P) = A exP[-Ag'(Pr/K)/RTl
(2)
where Ag' is the relative free energy as a function of PJp. Rewriting eq 2, we obtain Ag(P1/p) = -RT lnY(Pr/p)
+ A'
(3)
where A'is a constant. For our purpose, we need only the curvature of Ag'. Then, we put A' = 0 hereafter. We plot Ag' for the first shell as a function of P r / p in Figure 2. From this figure, we can clearly see that the curvature of the free energy curve increases very much with increase of z, as we e x p e ~ t e d . ' ~For - ~ ~the quantitative discussion, we define a relative curvature C by C = d / d , where d and w are the depth and width of the well, respectively, as shown in Figure 2. This curvature is calculated for many shells and is plotted in Figure 3 as a function of z2. It is seen that only the curvature of the 1st shell increases very rapidly and becomes almost in proportion to zz at large values
1010 The Journal of Physical Chemistry, Vol. 92, No. 5, 1988
2.0-
1.5-
”
I
I
I
/
/4 ,2
-5
2 0
1
2
3
22
Figure 3. Calculated values of the relative curvature Cas a function of z2 for some shells. The number in the figure denotes the shell number. The shell width is 1.1 A. The unit of C is arbitrary.
of z . The curvatures of the 2nd, 3rd, and 4th shells remain almost the same until z = 1 and increases a little when z > 1. For the other shells, the curvature is almost independent of z . When we calculate the free energy curve A i and the curvature by reducing the shell width to ca. 0.8 A, the 2nd shell exhibits some features of dielectric saturation for z = 1.0. This result suggests that the width of the dielectric saturation shell should be larger than 0.8 A. More detailed investigation has revealed that the shell width of the dielectric saturation is 1.1-1.2 A?2 The outside of this saturation shell constitutes the ordinary non-saturation region.
Discussion If we look at the saturation phenomenon from the molecular viewpoint based on the simulation data, seven molecules on the average are attached to the central reactant molecule and their positions are fluctuating in the radial direction with a width of about a half radius of the solvent molecule, which constitutes the shell width of the dielectric saturation, and an exchange of solvent molecules sometimes occurs between the one in the saturation shell and the one in the outside nonsaturation shell. The width of the saturation shell increases a little with increase of the valency of the molecular ion. The above results are in qualitative agreement with our previous phenomenological calculations based on the continuum model using the Lorentz field.23 In accordance with the dielectric saturation, the curvature of the free energy for the 1st shell as a function of solvent polarization increases very much with increase of the charge of the solute molecule and it becomes nearly proportional to z2 when z > 1. Indeed, we can prove that the curvature is proportional to z2 if the dielectric saturation is very strong, whereas it is a constant independent of z if no appreciable dielectric saturation occurs.22 The value of the relative curvature C is evaluated from Figure (22) Hatano, Y . ; Saito, M.; Kakitani, T.; Mataga, N., manuscript in preparation. (23) Kakitani, T.; Mataga, N. Chem. Phys. Lett. 1986, 124, 437.
Letters 3 to be 0.025 for z = 0 and 0.25 for z = 1. Therefore, the ratio of the force constant kJk, is 10, which leads to the @-value (defined by P = k,/(kc - k,)) of 0.1 1 . This P-value is in the range 0.1-0.3, which was extensively used in the analysis of the experimental data so far.l5-I9 As a check of the reliability of our calcuiations using a finite-size vessel, we have also conducted Monte Carlo simulations by eliminating solvent molecules fixed a t the surface of the vessel and/or by increasing the number N of the movable solvent molecules twice. In those calculations, the space-filling ratio is kept the same. The calculated results do not change from the original results for all the shells except a few shells near the surface of the vessel. This fact indicates that the dielectric saturation in molecular systems is quite a local phenomenon in the space just surrounding a charged molecule and it is little affected by the boundary condition. It will be instructive to discuss here the relation between the free energy curve Ag‘ and the total nonequilibrium free energy surface 3[P(r)] which is a functional of the polarization P(r) at distance r from the center. For practice, we have divided the continuum space into a large number of spherical shells. In each shell, the polarization fluctuates randomly and, accordingly, we can choose the polarization in each shell as an independent coordinate (or variable). Each free energy curve Ag‘ given by eq 3 should be plotted along the ith polarization coordinate. The total free energy 3 is actually given as a function in the multidimensional space of polarization coordinates,, The solvent equilibrium states for z = 1 and z = 0 will be represented by certain points in this space. This picture is the same as that given by Calef and W01ynes.~~The solvent polarization inside the saturation region is expected to be much more correlated to each other than in the outside regions. Therefore, it will be reasonable to treat the saturation region separately. The polarization coordinate in the saturation region is really the solvent coordinate of the c mode in our current treatment.l5-I9 Recently, Hwang and Warshel have carried out a molecular dynamics simulation for the ET reaction between benzene-like molecules in water, adopting the semiclassical trajectory m e t h d 2 * They have found that the solvent contribution to the activation barrier follows the Marcus’ relati~nship~,~ in the normal energy-gap region. It will be quite interesting to see by the similar simulation whether the Marcus’ relationship also holds true in the abnormal region where the effect of the dielectric saturation, if it takes place to any degree, will be demonstrated most clearly. On the other hand, it should be noted that the hydrogen-bond network of water molecules around the solute may play a significant role in disturbing the Occurrence of the dielectric saturation while such effect may be absent in the ordinary organic solvent such as acetonitrile.
Concluding Remarks The results of the present work using the spherical hard-core model for molecules show clearly that a thin dielectric saturation shell exists just around the charged molecule and that the curvature of the free energy curve due to the solvation increases very much with increase of the charge of solute, demonstrating the existence of the c mode. This c mode plays a crucial role in regulating the energy gap law of the ET, depending on the type of the reaction. The solvent vibration in the layers outside of the dielectric saturation shell constitutes an ordinary polaron mode called the s-mode in our previous treatment.Is-l9 Acknowledgment. T. K. was supported by a Grant-in-Aid (62580218) from the Japanese Ministry of Education, Science and Culture. The numerical calculations were done using FACOM VP-100 in Computation Center of Nagoya University and also using HITAC S-810/10 in the Institute for Molecular Science. (24) Calef, D. F.; Wolynes, P. G. J. Phys. Chem. 1983, 87, 3387. ( 2 5 ) Hwang, J. K.; Warshel, A. J . Am. Chem. SOC.1987, 109, 715.
