Reorganization energies for charge transfer reactions in binary mixtures of dipolar hard sphere solvents: A Monte Carlo study

J. Chem. Phys. 110, 473 (1999); https://doi.org/10.1063/1.478108 110, 473 © 1999 American Institute of Physics. Reorganization energies for charge transfer reactions in binary mixtures of dipolar hard sphere solvents: A Monte Carlo study Cite as: J. Chem. Phys. 110, 473 (1999); https://doi.org/10.1063/1.478108 Submitted: 07 July 1998 . Accepted: 25 September 1998 . Published Online: 21 December 1998 C. Denk, M. Morillo, F. Sánchez-Burgos, and Antonio Sánchez ARTICLES YOU MAY BE INTERESTED IN On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I The Journal of Chemical Physics 24, 966 (1956); https://doi.org/10.1063/1.1742723 On the Theory of Electron-Transfer Reactions. VI. Unified Treatment for Homogeneous and Electrode Reactions The Journal of Chemical Physics 43, 679 (1965); https://doi.org/10.1063/1.1696792 Electron–electron and electron-hole interactions in small semiconductor crystallites: The size dependence of the lowest excited electronic state The Journal of Chemical Physics 80, 4403 (1984); https://doi.org/10.1063/1.447218 Reorganization energies for charge transfer reactions in binary mixtures of dipolar hard sphere solvents: A Monte Carlo study C. Denka) and M. Morillo Universidad de Sevilla, Fı ´sica Teo ´rica, Apartado 1065, E-41080 Sevilla, Spain F. Sa ´nchez-Burgos and Antonio Sa ´nchez Universidad de Sevilla, Quı ´mica Fı ´sica, Facultad de Quı ´mica, C/ Profesor Garcı ´a Gonza ´lez s/n, E-41012 Sevilla, Spain ~Received 7 July 1998; accepted 25 September 1998! We study the behavior of the reorganization energy for simple charge transfer reactions in mixtures of dipolar hard sphere fluids by Monte Carlo simulation. The static dielectric constants of the solvents are also obtained from the simulation. They are used as input in the reorganization energy expressions provided by the Marcus theory and the mean spherical approximation. Thus, a comparison between the values obtained from the theoretical expressions and our simulation results is possible. The dependence of the reorganization energy with the mixture composition and the influence of preferential solvation effects is also discussed. © 1999 American Institute of Physics. @S0021-9606~99!51401-3# I. INTRODUCTION Electron transfer ~ET!reactions in polar solvents are important in many biological and chemical processes. Solvent fluctuations provide the transition state configurations for the solvent–solute system necessary for the ET. Marcus1–4 has given a physical picture for ET reactions, describing the multidimensional free energy surfaces of the product and reactant states as parabolic surfaces in terms of a suitable reaction coordinate. He related the activation free energy DG‡to the reaction free energy DGand the solvent reorganization energy l. Marcus also derived a simple expression for l, based on a macroscopic treatment of the solvent. The reorganization energy lis then given in terms of the static and optical dielectric constants of the solvent and a geometrical factor. Although the dielectric continuum model ~Marcus formula!provides an adequate qualitative description of the reorganization energy, its quantitative predictions are often times at variance with the experimental findings.5–7 Molecular descriptions of the solvent are, in principle, capable of overcoming some of the difficulties associated with a macroscopic treatment of the solvent. The mean spherical approximation ~MSA!treatment of the reorganization energy recognizes the solvent molecularity by describing the solvent molecules as hard spheres with point dipoles in their centers. Within the MSA, an expression for lhas been developed that includes the hard sphere radius of the solvent molecules.8,9 The need for a microscopic description is of particular interest in mixtures of polar solvents. For most mixtures one observes a nonideal ~nonlinear!behavior of the solvation energy of an ionic or dipolar solute with respect to the molar fractions of the species present in the solvent.10 This behavior is usually termed as preferential solvation. For ET reactions in mixtures of polar solvents, one might also expect an influence of preferential solvation effects on the solvent reorganization energy l. Both optical11 and thermal12,13 experimental data of ET reactions in mixtures of water and organic cosolvents show a behavior of lthat cannot be explained with a continuum model of the solvent. Analytical approximations for the microscopic description of the solvation energy in polar mixtures have been presented in the literature.10,14 In this paper we will investigate the behavior of the reorganization energy lin polar solvents by means of Monte Carlo ~MC!simulations. A quantitative prediction of experimental data12,13 would require very sophisticated simulation techniques. At this point we are rather interested in the general behavior of l. For this reason we have adopted a simple solvent model, with as few adjustable parameters as possible. Our solvent will be modeled either as a pure solvent or as a binary mixture. In both cases we will consider dipolar, nonpolarizable hard sphere molecules. The solute consists of two charged hard spheres. We will study charge separation and charge recombination processes for this system. As the solvent molecules are considered nonpolarizable, the optical dielectric constant will be e opt51 for all solvents. In order to compare the numerical results for lwith the theoretical predictions, we need to evaluate the static dielectric constant. Extensive simulations have been carried out in order to determine e 0for our different solvent models. The outline of the paper is as follows: In Sec. II we describe our solvent model and review the most important simulation details. In Sec. III we discuss the evaluation of the dielectric constants for the solvent models used in this work. In Sec. IV, the application of the MC technique to the calculation of free energy surfaces for charge transfer reactions is summarized. The results are presented and discussed in Sec. V. A comparison with the theoretical results as oba!Electronic mail: [email protected] JOURNAL OF CHEMICAL PHYSICS VOLUME 110, NUMBER 1 1 JANUARY 1999 4730021-9606/99/110(1)/473/11/$15.00 © 1999 American Institute of Physics tained for the continuum solvent model ~Marcus!and the microscopic description ~MSA!is made. Finally, Sec. VI contains some concluding remarks. II. MODEL AND MC SIMULATIONS Our solvent is modeled as a liquid of dipolar, nonpolarizable hard sphere molecules. Two solvent molecules iand k interact via a long ranged dipole–dipole interaction: Vi,k5 m i– m k rik 323~ m i–rik!~ m k–rik! rik 5rik>2rs~1! and a repulsive short ranged interaction ~the hard sphere potential, Vi,k5`for rik5 u rik u 5 u ri2rk u ,2rs). Here, rsis the radius of the hard spheres, m iis the dipole moment and riis the position of molecule i. In this work we study pure solvents and binary mixtures of dipolar hard sphere fluids. The radius of the solvent molecules is taken to be the same for all solvent components. A packing fraction of h 50.417 corresponding to a dense liquid was used in all simulations. The polarity of a pure solvent will be expressed in terms of the dimensionless parameter y54 p m 2 r 9kBT,~2! where r is the density of the liquid. A mixture of two solvents will be composed of two species, the less polar species L, characterized by its polarity yL, and another species H with a higher polarity yH. The composition of the mixtures will be described by the molar fractions of the more polar species fHthroughout this work. Periodic boundary conditions with the minimum image convention15,16 were applied to a cubic simulation box of side length L. This method requires the use of a cutoff for the long ranged dipolar interactions. In order to account for the contributions beyond the cutoff radius rc<L/2 we have adopted the generalized reaction field method.17 We consider the moving boundary dielectric implementation, i.e., the subsystem inside the cutoff sphere around each molecule is thought to be immersed in a continuum dielectric characterized by a macroscopic static dielectric constant e RF . A molecule interacts directly with all molecules inside its cutoff sphere via Eq. ~1!and with the reaction field. The reaction field at the center of molecule idue to the molecules inside the cutoff sphere surrounding it is proportional to the total dipole moment inside its cutoff sphere16 m i–Ei r52~ e RF21! 2 e RF11 1 rc 3 m i•( rik<rc m k5 am i•( rik<rc m k,~3! where we have defined a screening constant a . Summing up all terms contributing to the total potential energy we obtain Etot51 2( kÞi,rik<rc Vi,k dd 21 2 a ( i m i 2~4! with an effective interaction energy Vi,k dd 5 m i– m k~12 a rik 3! rik 323~ m i–rik!~ m k–rik! rik 5.~5! The last term in Eq. ~4!describes the constant selfinteraction of the molecules with their own reaction fields. The form of the effective interaction in Eq. ~5!is very convenient for computational purposes. The inclusion of the reaction field only requires two extra floating point operations for each evaluation of the interaction energy, which renders this method computationally feasible, even for relatively large systems. It has been shown that the computationally more expensive Ewald summation technique yields equivalent results for the dielectric constant when applied to similar systems.18,19 Analogous conclusions have been drawn for free energy calculations of ionic hydration.20 System configurations were generated in the following manner: starting from a given configuration, a molecule was selected randomly. It was displaced from its initial position in a cube of side length Dxwith a uniform distribution. Now, a rotational axis (x,yor z) was selected at random and the dipole orientation of the molecule was rotated by a uniformly distributed angle in the range 2D f < f <D f around this axis. The parameters Dxand D f were adjusted to reach an acceptance ratio of approximately 30%. The implementation of the simulation is along the lines of the standard Monte Carlo techniques.15,16 III. DIELECTRIC CONSTANT The application of the statistical mechanical theory of the dielectric constant21,22 to finite size simulation systems with boundary conditions requires some modifications. This issue has been lucidly addressed by Neumann.18 From his analysis it follows that the static dielectric constant e 0in the reaction field ~RF!geometry is given by e 052 e RF~11 z !11 112 e RF2 z ,~6! with z 54 p 3 b ^ M2 & L353ygK.~7! Here ^ M2 & is the configurational average of the squared total dipole moment of the system and b 51/kBT. We have also indicated the relation of z with the polarity yand the Kirkwood gfactor gK5 ^ M2 & /N m 2. The dielectric constant e RF that characterizes the reaction field was determined selfconsistently in the simulations, so that e RF' e 0for all solvents. The solvent molecules were initially prepared in an fcc structure and the system was then allowed to relax during 53106MC configurations. Mean values were obtained from NMC553108subsequent MC configurations, except where otherwise stated. In order to obtain an estimate of the statistical errors, we calculated mean values of z over blocks of 106configurations. We carried out various tests to check the convergence of the values of the dielectric constant obtained in our simulations. In Fig. 1 we show the running averages of e 0for two solvents with y52.18 and y52.90 for two different values of the maximum rotational angle D f . In the case of the solvent with the higher polarity, the values of e 0converge very slowly. Even for NMC553108configurations we 474 J. Chem. Phys., Vol. 110, No. 1, 1 January 1999 Denk et al. still find a deviation of D e 051.6 between the two simulations ~theoretically, the simulation results should be independent of D f ). The statistical errors estimated using the block averages were D e 050.4 for simulation ~a!and D e 050.5 for simulation ~b!. This indicates that, even with 53108MC configurations, the phase space of this highly polar system has not been sufficiently explored to obtain precise values of e 0. The deviation of the two curves will be used as an estimate for the error of the dielectric constant ~approximately 3%!for solvents with high polarity. Solvents with lower molecular dipole moments show a much faster convergence of e 0. For solvents with a dielectric constant e 0<30 we found 23108MC configurations to be sufficient to determine e 0 with a relative error of approximately 3%. We have also studied the dependency of e 0on the system size. Simulation ~a!of Fig. 1 was repeated for a system with N5864 particles; the resulting value of e 0is shown in Table I @simulation ~e!#. The results indicate that e 0is independent of the size of the system for N>256 within the error limits. As already mentioned, the parameter e RF was adjusted self-consistently by repeating each simulation, using the result e 0of a simulation as an input parameter e RF in the next simulation. In all cases, it was sufficient to repeat each simulation only once, as e 0depends very weakly on e RF .In simulation ~f!~see Table I!we have repeated simulation ~b! with e RF558. The results are practically identical. The dielectric constant e 0was then determined for a pure solvent over a wide range of values of the molecular polarity y. A system with N5256 particles was used in these simulations and the mean values of e 0were obtained after 53108MC configurations. For the highest polarity, y 53.0, we increased the number of MC configurations to NMC5109for the reasons mentioned above. The results are shown in Fig. 2. The theory of liquids provides various approximations for the structure of the dipolar hard sphere fluid.23 Of particular interest is the so-called MSA, as it provides analytical expressions for the correlation functions and the dielectric constant.24 We find it instructive to compare our simulation results with the MSA theoretical predictions. The pair distribution function can be expanded as23 h~1,2!5hS~R!1hD~R!D~1,2!1hD~R!D~1,2!,~8! where D~1,2!is the cosine of the angle formed by the dipole orientations of two molecules and 2D(1,2) m 2/R3is the dipole–dipole interaction as defined in Eq. ~1!. The MSA provides the functions hS(R),hD(R) and hD(R) in terms of the radial distribution function ~rdf!of the Percus–Yevick ~PY!solution for hard spheres at different densities. In Fig. 3 we show the radial distribution function gS(R)5hS(R)11 for a highly polar solvent (y53.00). Solvents with a lower polarity have a very similar rdf with a slightly lower main peak. The MSA result for gS(R)~also shown in Fig. 3!is just the PY rdf for hard spheres at density r and does not depend on the molecular polarity. Although the agreement is globally good, there are deviations between the MSA and the simulations. In the region close to contact, the MSA subestimates the value of gS(R), which is of primary importance for many thermodynamic properties of the liquid. It also predicts a somewhat slower decay from the peak value to the first minimum when compared with the simulations. The dielectric properties of the solvent are related to hD(R), which describes the angular correlation of two molecules at a given distance R. The dielectric constant e 0is FIG. 1. Running averages of e 0for ~a!,~b!y52.90 and ~c!,~d!y52.18. The systems were simulated with two different values of the maximum rotational angle D f 5 p /2 ~a!and ~c!and D f 51.0 ~b!and ~d!. FIG. 2. The static dielectric constant e 0for a pure solvent. The error bars of the simulation results indicate the estimated relative error of 3% ~see main text!. The solid line is a fourth degree interpolation polynomial. The dotted line corresponds to the theoretical result as obtained from the MSA. TABLE I. Dielectric constants obtained for solvents with y52.90 ~a!,~b!,~e!,~f!and y52.18 ~c!,~d! e 0NMC ND fe RF ~a!59.68 53108256 p /2 70 ~b!58.05 53108256 1.0 70 ~c!28.29 53108256 p /2 30 ~d!29.05 53108256 1.0 30 ~e!58.41 23108864 p /2 70 ~f!58.78 53108256 1.0 58 475J. Chem. Phys., Vol. 110, No. 1, 1 January 1999 Denk et al. determined by the polarity yand the Kirkwood gfactor gK @see Eq. ~6!#. The Kirkwood gfactor gKcan be obtained from hD(R) by integration gK5114 pr 3 E 0 `hD~R!R2dR.~9! In Fig. 4 we show hD(R) for a highly polar solvent (y 53.0) as obtained from our simulations and the corresponding MSA result. The MSA does not provide a good approximation to hD(R) as it underestimates the angular correlation. This results in an underestimation of gKand thus e 0. This deviation is more pronounced for solvents with a high polarity ~see Fig. 2!. But, even for solvents with the lowest polarity under consideration in this work, the MSA result for hD(R) does not agree with the simulation data. The dielectric constant as obtained by the MSA gives a good approximation to e 0only for solvents with gK'1. These limitations of the MSA are well known and more sophisticated theories based on the hypernetted chain ~HNC!approximation provide much better angular correlation functions and dielectric constants for dipolar hard sphere fluids. Comparisons of the MSA and other theories can be found in the literature.23,25–33 Let us now consider the case of binary solvent mixtures. We study two different types of mixtures: ~A!mixtures with two components of similar polarities (yH53.0 and yL 52.18), and ~B!mixtures whose components are of rather different polarity (yH53.0 and yL50.75). For each type of mixture, the dielectric properties will only depend on its molar composition ~defined by the molar fraction of species H, fH), if the temperature is held fixed. This gives us the possibility of studying two quite different mixtures, spanning a wide range of dielectric constants. The simulation procedure is the same as for pure solvents, and the dielectric constant is evaluated by using Eqs. ~6!and ~7!. In Fig. 5 we show the dielectric constants as obtained for the simulated compositions of mixtures of type ~A!and ~B!. The values of e 0for the pure solvents (fH50 and fH51) were taken from the simulations described above ~see Fig. 2!. For both mixtures we observe an almost quadratic dependence of the static dielectric constant on the molar fraction fH. This behavior can be understood by noticing that in the case of a binary mixture, ^ M2 & /Nis to a good approximation a quadratic function of the molar fraction fH. For e 0 ' e RF , it follows from Eq. ~6!that the dielectric constant is essentially a linear function of z , which in turns is proportional to ^ M2 & /Nas can be seen in Eq. ~7!. IV. REORGANIZATION ENERGIES A. Calculation of the reorganization energy from molecular simulations In Sec. III we have analyzed the dielectric behavior of dipolar hard sphere solvents. We will now study the situation when a solute is immersed in the solvent. Our aim here is to study the energetics of thermal charge transfer reactions between two solute molecules in the presence of a polar solvent. We have adopted a simple model for the solute: The solute consists of two hard sphere molecules ~donor and acceptor!with given radii rdand raseparated by a fixed disFIG. 3. The radial distribution function gS(R)(R5r/2rs) for a solvent with y53.00 ~filled circles!. The dotted line represents gS(R) as obtained from the MSA. FIG. 4. The angular correlation function hD(R) for a highly polar solvent (y53.0, filled circles!. The dotted line corresponds to hD(R) as provided by the MSA. FIG. 5. Static dielectric constant e 0for binary mixtures of dipolar hard spheres with yL52.18 @mixture ~A!, circles#and yL50.75 @mixture ~B!, squares#in terms of the molar fraction fH. In both cases yH53.0. The solid lines represent quadratic interpolation polynomials. 476 J. Chem. Phys., Vol. 110, No. 1, 1 January 1999 Denk et al. tance d. The donor and acceptor in the reactant state carry a charge qdand qa, respectively; a negative net charge transfers from the donor to the acceptor. We will consider the reaction coordinate of this process in terms of a charging parameter j , so that qd~ j !5qd1 j e,~10! qa~ j !5qa2 j e.~11! Here, eis the elementary charge ~that we will take as positive!. The reactant state of the solute is given by j 50, while the product state is obtained by setting j 51. In this work we have studied two typical transfer reactions: ~i!a charge recombination ~CR!process qd52eand qa51eand ~ii!the inverse process, a charge separation ~CS!process (qd5qa 50). The rates of these processes are governed by an activation free energy law. Marcus1–4 obtained for the nonadiabatic electron transfer rate ket5 k A 4 p l/ b exp~2 b DG‡!,~12! where the activation free energy DG‡is defined as DG‡5~DG1l!2 4l.~13! Here, k is a matrix element describing the electronic coupling between reactant and product state, and DGand l denote the free-energy change of the reaction and the reorganization energy, respectively. Equation ~12!was derived under the assumption of a classical solvent that responds linearly to a redistribution of charges. The procedure used for the evaluation of reorganization energies in charge transfer reactions from simulations is well documented.34–36 Here we will briefly indicate the main points. Let H j denote the solvent–solute interaction energy for a solute state characterized by the charging parameter j and a certain fixed solvent configuration. The energy gap DV5H12H0describes the energy difference between products and reactants for a given configuration of the solvent. An electron transfer will take place only for solvent configurations that fulfill the condition DV52Ei~Frank–Condon principle!, where Eiis the intrinsic energy difference between the gas phase electronic structures of the initial and final state of the solute. Let us define the random variable D5DV, with a probability law given by p j (D) 5 ^ d (D2DV) & j , where the angular brackets indicate an equilibrium average taken with a canonical distribution describing a system at temperature Tand Hamiltonian H j . The rate of the ET process is proportional to the probability density p0(2Ei) of the energy gap DVhaving a value 2Ei. This probability distribution corresponds to a solvent in thermodynamical equilibrium with the reactant state of the solute, j 50. The main problem in simulations lies in the construction of this probability density as the entire phase space of the solvent degrees of freedom has to be explored. This is an impractical task, and one resorts to a free energy perturbation method,34 which is based on the following. For a given value of the fractional charge parameter j ,DVis sampled around a value D j with a distribution that is approximated very well by a Gaussian p j ~D!51 A 2 ps j 2e2~D2D j !2/2 s j 2.~14! Various states j of the solute are simulated, each simulation providing a probability distribution p j (D) around D j . These distributions can be pieced together34 to yield the distribution p0(D) over a wide range of values D. Zhou and Szabo35 have proposed a simplified method that permits obtaining p0(D) by simulating only two states of the solute: j 50 and j 51. In their work, the reorganization energy is given by l5D02DG10 ,~15! with the free-energy difference between state j and the reactant state of the solute DG j 05 E 0 j d j 8D j 8.~16! The mean value of the energy gap D j for a certain value of the charging parameter j can be obtained by a cubic interpolation polynomial in j , where all coefficients are determined by the mean values D0and D1and the dispersions b s 0 2and b s 1 2in the reactant and product state, respectively. Inserting the polynomial expression of Zhou and Szabo in Eqs. ~15!and ~16!leads to lCR51 2~D02D1!11 12 ~ b s 0 22 b s 1 2!,~17! lCS51 2~D02D1!21 12 ~ b s 0 22 b s 1 2!.~18! As pointed out by Zhou and Szabo, the accuracy of the interpolation polynomial approximation for D j can be checked by carrying out an additional simulation for j 51/2. The value of the energy gap D1/2 obtained is then compared to the value provided by the interpolation formula. We found a deviation of less than 1% for all cases considered here. From the interpolation polynomial, the probability density p0(D) can also be constructed.35 B. Simulation details The inclusion of the charged solute particles in the simulation introduces a difficulty when dealing with systems of finite size: The orientation of the dipole moments of the solvent molecules is anisotropic and the periodic replication of the central simulation cell does not represent a physical picture of the solvent. This issue has been widely discussed in the literature37 and various methods have been developed to avoid periodic boundary conditions ~pbc!for such systems. Spherical simulation cells in conjunction with a suitable method to avoid self-polarization on the surface of the cell are usually employed.38 We have considered simulations with periodic boundary conditions and within a spherical geometry. In the pbc simulations, the charge–dipole interactions can be described via an effective interaction potential that includes the reaction field in a similar manner to Eq. ~5!: 477J. Chem. Phys., Vol. 110, No. 1, 1 January 1999 Denk et al. Vi,j qd5qi S 2 rij 2 a D m j–rij,~19! where a is the screening constant as defined in Eq. ~3!. The solute charges interact with all dipoles inside their corresponding cutoff spheres. The cutoff radius for the charge– dipole interactions was set to be the same as for the dipole– dipole interactions. In the spherical model, the simulation cell is a spherical vessel with the solute situated at the center of the vessel. Each molecule interacts with all molecules inside the vessel and electrostatic interactions are taken fully into account @set a 50 in Eqs. ~5!and ~19!#. The simulation sphere is divided into two regions: an inner sphere, where the molecules are allowed to move and rotate and an outer shell, where the molecules are kept fixed during the simulation. The molecules constituting the outer shell are situated initially in a fcc structure with randomly distributed dipole moments. The fixed dipoles in the outer shell prevent an unphysical polarization of the dipoles close to the surface of the vessel. A correction to the energy gap D0has to be applied in the case of the spherical model. Contributions from the solvent outside the sphere can be partially taken into account by considering the sphere to be immersed in a continuous dielectric medium with dielectric constant e out . The influence of this dielectric medium is not taken into account during the simulation ~as it would require the full solution of the Laplace equation inside the sphere for each configuration!, but its effect on D0may be estimated by electrostatic considerations. For a CR process, the contributions to D0from outside the simulation cell with radius Rare D0 out52~ e out21! 2 e out11 ~ed!2 R3,~20! where dis the distance between the solute charges. We studied the dependency of D0on the size of the system for both models. A highly polar solvent (y52.90) was used to solvate two solute ions with equal radii rd5ra 5a51.44Å, charges qd52e,qa51eand separated by a distance of d56 Å. The solvent radius was set to be the same as the solute radius and a packing fraction of h 50.417 was used. In the case of the pbc geometry we simulated systems with N5246 490 854 solvent molecules. The cutoff radius for the dipole–dipole and charge–dipole interactions was set to rc5min(L/2,8rs) for each system, where Lis the sidelength of the cubic simulation cell. For the system with N 5854 solvent molecules we conducted a simulation with the full cutoff rc5L/2 in order to check the influence of the reduced cutoff on D0. In the spherical geometry we simulated the same system with a total of Ntot5673 solvent molecules. In this case, the radius of the inner sphere ~where the molecules are allowed to move!was varied. The results are combined in Table II. For small system sizes, the two geometries produce quite different values for the energy gap D0. The pbc simulations overestimate D0, whereas the spherical model results in too small values for the gap. As the system size increases, both models yield similar results. As the size of the outer shell in the spherical model decreases ~last value in Table II!, the polarization of the cell surface results in lower values for D0. The reduced cutoff rc58rs511.52Å for the pbc system with N5854 solvent molecules gives practically the same value for D0as the system with rc5L/2514.75Å. Although the spherical geometry permits the usage of smaller system sizes, we have opted for the pbc geometry with N5854 solvent molecules and a reduced cutoff rc 58rsfor our simulations. Our choice is based on the fact that the spherical model was not able to reproduce the dielectric constant of the solvent. Thus the usage of a single geometry for the determination of both the dielectric constant and the reorganization energy is clearly preferable. The pbc geometry also allows us to calculate the solvent radial distribution functions. With our choice of the system size and the cutoff radius, a molecule close to the simulation cell boundary only interacts with a small fraction of the periodic images of the solvent molecules in the first solvation shell. When applying the conventional Monte Carlo method for creating system configurations ~translation and rotation of a molecule!in the case of mixtures, we noticed a very slow relaxation of the system. If a charged solute is present, the initial fcc structure of the solvent relaxes rapidly to a situation where the solute is solvated by the molecules that are initially in the vicinity of the solute. Thus, the molar composition of the first solvation shells depends on the initial positions of the two species within the fcc grid. The high density of solvent molecules close to the solute now renders a restructure of the solvation shell extremely difficult. If the two components of the mixture are equally sized, one can adopt a much more efficient method for creating configurations: The conventional ‘‘move’’ is alternated with a ‘‘swap’’ of molecules: two molecules of different species are selected randomly. Now, the identities of the molecules are interchanged, i.e., the dipole moment of the two molecules is changed. We only change the absolute values of the dipole moments; the dipole orientations are not altered. Each new configuration is generated either by a move or by a swap and the probability for a swap was set to 10%. V. RESULTS AND DISCUSSION The reorganization energies lCR and lCS for a charge recombination and a charge separation process have been TABLE II. Values of D0for pbc and the spherical model. The solvent and solute parameters are specified in the main text. For the spherical model N refers to the number of molecules that are allowed to move. In this geometry the total number of molecules was held fixed at Ntot5673. ND0~kcal/mol!rcGeometry 246 226.9 rc5L/2 pbc 490 223.1 rc58rspbc 854 221.9 rc58rspbc 854 222.1 rc5L/2 pbc 239 217.9 ¯spherical 371 222.4 ¯spherical 521 223.5 ¯spherical 593 218.6 ¯spherical 478 J. Chem. Phys., Vol. 110, No. 1, 1 January 1999 Denk et al. calculated for various solvents using the methods described above. The solute consists of two hard spheres with radius rd5ra5a51.44Å and distance d56Å. A. Pure solvents We first study the case of ET reactions in pure solvents. We have considered solvents with polarities yin the range of 0.75<y<3.0. The solvent radius was chosen to be the same as the solute radius. A packing fraction of h 50.417 and a temperature of T5300K was used in all simulations. We will compare our simulation data with the reorganization energy expressions provided by the Marcus theory and the MSA. In the Marcus continuum description, the reorganization energy for a charge transfer of an elementary charge between two equally sized solute molecules of radius aat a distance din a nonpolarizable solvent with static dielectric constant e 0is given by l5e2 S 121 e 0 DS 1 a21 d D .~21! The MSA theory takes into account the molecular aspect of the solvent. In the case of infinitely separated charge centers (d→`), the reorganization energy is just the sum of the solvation free energies of the two ions8,9 l5e2 S 121 e 0 D 1 a~11 d !.~22! Here, d is a correction to the solvated ion radius, which to a very good approximation is given by39 d 53rs/a~1081/3 e 0 1/622!21.~23! If the distance between the ions is not infinite, the inclusion of the ion–ion interaction in Eq. ~22!within the framework of the MSA would require the knowledge of the mean ion– ion potential. We will approximate this term by the screened interaction in a continuum dielectric medium, as it has been done by others.36 The inclusion of this term in Eq. ~22!results in lMSA5e2 S 121 e 0 DS 1 a~11 d !21 d D .~24! We should keep in mind that this approximate lMSA does not emerge from a purely microscopic picture of the solvent, as its influence on the ion–ion interaction is taken into account by a continuum description. This is expected to be a good approximation, as long as the distance of the ions is sufficiently large. When evaluating Eqs. ~24!and ~23!, we have to provide the solute radius a, the ion distance d, the solvent radius rs and the static dielectric constant of the solvent e 0. The latter may be obtained directly from the MSA, using the following relations:23 3y5~114 j !2 ~122 j !42~122 j !2 ~11 j !4, ~25! e 05~114 j !2~11 j !4 ~122 j !6. For comparison with experimental data it is often more convenient to use the experimental value of e 0in Eqs. ~24!and ~23!. We will thus compare three theoretical expressions with our simulation data: ~a!the Marcus expression Eq. ~21!,~b! the ‘‘consistent’’ MSA result as given by Eq. ~24!using the solvent polarity yin order to determine e 0from the MSA and ~c!the ‘‘experimental’’ MSA result lMSA ex , using the dielectric constants as obtained from the simulations in Eq. ~24!. The results are shown in Fig. 6, where we represent the reorganization energy versus the Pekar factor 121/ e 0. In this representation, the Marcus result @Eq. ~21!# is a straight line. We find that lCS.lCR for all simulated pure solvents. Similar results have been obtained by other authors for comparable systems.36 This is due to the combination of two effects:34,40 First, the free energy surfaces are not strictly parabolic, as it would be required by linear response theory. Second, even if the deviations from parabolic surfaces are small, the curvature of the parabolas may be different in the reactant and product states. Both features give rise to a different reorganization energy for the charge separation and charge recombination processes. The Marcus expression for lresults in an overestimation of the reorganization energy. Equation ~21!gives l(y 53.0)5172.6 kcal/mole for the solvent with the highest polarity. In the plot we have scaled the corresponding curve to coincide with lMSA ex at this polarity. This is equivalent to using an effective ion radius of a52.0Å in Eq. ~21!. When the geometric factor is adjusted in this way, the continuum description still does not describe well the overall behavior FIG. 6. Reorganization energies for a pure solvent. The squares represent lCS, the circles correspond to lCR as obtained from the simulations. The dashed line is the result as obtained by the continuum description ~Marcus! scaled to coincide with lMSA ex for the highest polarity. The dotted line corresponds to the theoretical MSA result for the simulated system. The MSA result lMSA ex , using e 0as obtained from the simulations as an input parameter, is represented by the solid line. 479J. Chem. Phys., Vol. 110, No. 1, 1 January 1999 Denk et al. of the reorganization energy. This feature is expected. The expression for lin the continuum description, Eq. ~21!, indicates that the reorganization energy depends only on e 0for fixed radii and separation distance and this dependence is very weak. For the MSA results, we find a much better agreement between the theoretical expressions and the simulation data. The theoretical MSA result ~b!gives good agreement for small polarities, but underestimates the reorganization energy for high polarities. This is not surprising, as it is well known that the MSA underestimates the dielectric constant for high solvent polarities. This results in an underestimation of the Pekar factor 121/ e 0on one hand and of the geometric factor 1/a(11 d ) on the other. The MSA expression lMSA ex with e 0taken from the simulation data yields reorganization energies which are close to the simulation results for the charge recombination process. By construction, the MSA does not distinguish between the charge recombination and charge separation processes, which as shown by the simulation results, have different reorganization energies. Thus, the good agreement of lMSA ex with the results for the charge recombination process in the case presented is probably fortuitous. However, this expression describes the overall shape of both curves much better than the Marcus approximation. The theoretical expression for lwithin the framework of MSA includes the geometric factor 1/a(11 d ) which depends on e 0via Eq. ~23!. The ability of lMSA ex to reproduce the trend of the simulation data for both processes suggests that the inclusion of this geometric factor is an important correction to the Marcus expression for the reorganization energy, Eq. ~21!. The geometric factor, as given by the MSA, can be interpreted as an effective solute radius a(11 d ) and is a consequence of the solvent molecularity in the MSA picture. This emphasizes the need for a molecular description of the solvent, even in the case of pure solvents. From these results one may be tempted to conclude that the MSA is able to predict the solvent structure around the solutes. We are currently examining this question for single ions in solution. Preliminary results indicate that the polarization density of the solvent around an ion is not very well represented by the MSA expression. However, the free energy of solvation, given as an integral over the polarization density, is predicted quite well by the MSA expression. B. Mixtures We will now proceed to the case of binary mixtures of polar solvents. We have calculated the reorganization energy for the two mixtures studied in the previous section @mixture ~A!with similar polarities of the two components: yH53.0 and yL52.18 and mixture ~B!with rather different polarities of the two components: yH53.0 and yL50.75#. The results for mixture ~A!are shown in Fig. 7 and for mixture ~B!in Fig. 8. When the solute is not preferentially solvated by one of the two components of the mixture, one would expect a linear behavior of the solvation energy, as the molar fraction of the components is varied.10 For nonpolarizable solvents, the theoretical expressions for the reorganization energy reduce to a solvation energy in the limit d→`. Thus, a linear ~or close to linear!dependence of lon the molar fraction fH would indicate that no preferential solvation is present in the process. For both mixtures we observe deviations from a linear behavior in the reorganization energy. In mixture ~A! we observe an almost linear behavior for lCS, while all values of lCR for the mixtures lie above a straight line connecting the results for the pure solvents. For mixture ~B!this behavior is much more pronounced. When adding a small molar fraction of the component with the higher polarity, the reorganization energy shows a drastic increase, and very quickly the reorganization energy lreaches a value close to that of the pure solvent (fH51). As it was already observed for mixture ~A!, the reorganization energy for the charge recombination process lCR deviates more from the linear FIG. 7. Reorganization energies for mixture ~A!(yH53.0 and yL52.18). The squares represent lCS, the circles correspond to lCR. The linear dependence, as expected in the limit d→`for ideal solvation, is represented as a dashed line for each case. The xsymbols represent the MSA results lMSA ex . FIG. 8. Reorganization energies for mixture ~B!(yH53.0 and yL50.75). The squares represent lCR, the circles correspond to lCS. The solid lines connecting the data points have been obtained by spline interpolation. The linear dependence, as expected in the limit d→`for ideal solvation, is represented as a dashed line for each case. The xsymbols represent the MSA results lMSA ex . 480 J. Chem. Phys., Vol. 110, No. 1, 1 January 1999 Denk et al.