Solvation in binary mixtures of dipolar hard sphere solvents: Theory and simulations

J. Chem. Phys. 113, 2360 (2000); https://doi.org/10.1063/1.482051 113, 2360 © 2000 American Institute of Physics. Solvation in binary mixtures of dipolar hard sphere solvents: Theory and simulations Cite as: J. Chem. Phys. 113, 2360 (2000); https://doi.org/10.1063/1.482051 Submitted: 16 February 2000 . Accepted: 10 May 2000 . Published Online: 28 July 2000 M. Morillo, Claus Denk, Francisco Sánchez-Burgos, and Antonio Sánchez ARTICLES YOU MAY BE INTERESTED IN Invariant Expansion for Two-Body Correlations: Thermodynamic Functions, Scattering, and the Ornstein—Zernike Equation The Journal of Chemical Physics 56, 303 (1972); https://doi.org/10.1063/1.1676864 Invariant Expansion. II. The Ornstein-Zernike Equation for Nonspherical Molecules and an Extended Solution to the Mean Spherical Model The Journal of Chemical Physics 57, 1862 (1972); https://doi.org/10.1063/1.1678503 Invariant expansion III: The general solution of the mean spherical model for neutral spheres with electostatic interactions The Journal of Chemical Physics 58, 3295 (1973); https://doi.org/10.1063/1.1679655 Solvation in binary mixtures of dipolar hard sphere solvents: Theory and simulations M. Morillo and Claus Denk Universidad de Sevilla, Fı ´sica Teo ´rica, Apartado 1065, E-41080 Sevilla, Spain Francisco Sa ´nchez-Burgos and Antonio Sa ´nchez Universidad de Sevilla, Quı ´mica Fı ´sica, Facultad de Quı ´mica, C/ Professor Garcı ´a Gonza ´lez s/n, E-41012 Sevilla, Spain 共Received 16 February 2000; accepted 10 May 2000兲 The structure of mixtures of dipolar hard sphere fluids with components of equal size but different dipole moments around a single ion is studied. The solvation energy and the polarization around the ion is obtained in the framework of the mean spherical approximation 共MSA兲. Our theoretical results and the results of other workers are compared with simulation data obtained from Monte Carlo simulations. An interpretation of the meaning of preferential solvation is given in terms of the contrasting behaviors of partial polarization in the bulk and near the ion. © 2000 American Institute of Physics. 关S0021-9606共00兲50830-7兴 I. INTRODUCTION The solvation of ions by polar solvents has been the subject of numerous theoretical and experimental studies. One of the fundamental quantities is the free energy of solvation. Its simplest description relies on the consideration of the solvent as a macroscopic dielectric continuum that becomes polarized by the presence of the ion. The interaction energy between the solvent polarization and the field created by the ion is then characterized by the dielectric constant of the solvent. Although the macroscopic approximation captures the basic ingredients of the free energy of solvation, their quantitative predictions are often at variance with the experimental findings, in particular for solute ions whose sizes are not very large compared with the size of the solvent molecules. The limitations of the dielectric continuum treatment of the Born free energy of solvation has led to the introduction of other alternatives that contemplate the molecular description of the solute–solvent system. Starting from the integral equations of liquid theory, and under suitable approximations, analytical expressions for the solvation free energy have been found. In particular, for a model of a hard sphere solute ion in a bath of solvent molecules formed by hard spheres with point dipoles, Chan et al. derived a formula for the Born free energy of solvation within the mean spherical approximation 共MSA兲.1Their expression amounts to replacing the solvent hard sphere diameter in the classical Born expression by an effective one, which depends on the dipole number density, the dipole moment, and the dielectric constant of the pure solvent. Their theoretical treatment gives support to empirical expressions widely used to fit experimental data. A molecular theory of solvation based on the density functional theory has been put forward by Chandra and Bagchi,2leading to an expression for the free energy of solvation which has to be evaluated numerically. Patey and co-workers3have used the linearized hypernetted chain closure 共LHNC兲to solve the Ornstein–Zernike equation for the same model of hard sphere dipoles embedding a hard sphere point charge. The resulting expressions are more complicated than the corresponding ones found in the MSA, their solutions requiring a nontrivial numerical treatment. It is found that the solvation energy obtained within the LHNC agrees well with that obtained from the MSA, except for solute ions whose diameters are much larger than that of a solvent molecule. Other theoretical approaches are based on the use of perturbation theory. In recent work, Pade ´approximant techniques have been exploited to truncate the perturbation expansion of the solvation chemical potential of a dipole in dipolar liquids.4Extensions of the dielectric continuum theory that include dielectric saturation and electrostriction effects have been developed in Ref. 5. Numerical simulations have also been profusely used in the study of solvation problems.6–12 General considerations about the calculation of free energies of solvation can be found in Ref. 13. Although in recent years most of the interest has focused on the analysis of the dynamics of solvation,14–18 some work has been devoted to the study of structural and thermodynamical properties. Monte Carlo simulations were used in Ref. 6 to analyze the dependence of the solvation energy in a polar liquid with the solute charge and the influence of dielectric saturation effects. In Ref. 9, the dependence on the ionic charge of hydration free energies of ions is studied. A detailed analysis of the thermodynamics of ion solvation in dipolar fluids using Monte Carlo simulations and the mean reaction field method is presented in Ref. 12. Recently, simulations of a dipole in a bath of polar hard spheres have been presented in Ref. 4. In this paper we address the problem of calculating the solvation free energy of an ion in a mixture of polar fluids. Mixtures of polar solvents are very convenient from an experimental point of view as the polarity of the solvent can be easily controlled by varying the composition of the mixture. Nonetheless, the analysis of solvation of ions in mixtures has received less attention than in pure solvents.2,19 In a previous JOURNAL OF CHEMICAL PHYSICS VOLUME 113, NUMBER 6 8 AUGUST 2000 23600021-9606/2000/113(6)/2360/9/$17.00 © 2000 American Institute of Physics study,20 we analyzed the dependence of the reorganization energy of electron transfer reactions in polar mixtures with the mixture composition, using Monte Carlo techniques. The simulations clearly indicate that there exists an excess reorganization energy: as a small amount of the solvent with higher polarity is added to the mixture, the reorganization energy increases drastically. There, the inability of theoretical expressions to explain this phenomena was noticed. The calculation of the solute–solvent radial distribution function seemed to indicate that the microscopic origin of this phenomenon is the preferential solvation of the ions by the more polar species. The theoretical study of the influence of preferential solvation on the reorganization energy is complex, as one has to deal with two charge centers. Here, we consider the somewhat simpler problem of a single ion in a solution at infinite dilution and carry out a theoretical and simulation analysis of a model system. We will see that a key ingredient to understand the behavior of the free solvation energy of an ion in a mixture of polar solvents is the preferential solvation of the ion by the component with the higher polarity. The existence of preferential solvation has been related in the literature to deviations of the free solvation energy from an ideal behavior, characterized by a linear dependence of the free energy with the mixture composition.2We believe that deviation from that assumed linearity is not a signature of preferential solvation. Rather, we focus on the behavior of the total and partial polarization densities of the two components around the ion induced by the ion fields. They contain the key solute– solvent correlations and, based on their knowledge, one can give a microscopic interpretation of what is normally termed preferential solvation. Namely, the composition of the mixture near the ion is substantially different from the nominal composition of the liquid in the bulk, away from the ion. This feature can be characterized by the excess local molar fraction parameter 共Marcus兲.21,22 The component with the higher polarity contributes to the total polarization density of the first two solvation shells around the ion with a much greater value than its bulk molar fraction would indicate. This local structuration of the mixture in the first few solvation shells gives rise to preferential solvation, and as the authors of Ref. 2 point out, its existence is not necessarily related to specific details of the molecular interactions such as hydrogen bonds. Furthermore, we notice that continuum dielectric theories for the free solvation energy containing the dielectric constant as input are not valid, in general. The reason is that the dielectric constant incorporates just the solvent–solvent correlations in the absence of solute, and these correlations are not the relevant ones in relation with the solvation free energy. The structure of the paper is as follows. In Sec. II, we present analytical results for a model of a hard sphere solute ion embedded in a mixture of hard sphere dipolar solvents of equal radii and different polarities. The MSA theory for the case of single component solvent has been thoroughly developed by Chan et al.1We extend their ideas to the case of a binary solvent mixture, obtaining analytical expressions for the partial and total solvent polarizations and the free energy of solvation. We also briefly review the theory of Chandra and Bagchi2based on a linearization of the free-energy functional of an inhomogeneous mixture in the presence of an ion field. In Sec. III we describe the methodology used in our Monte Carlo simulations. Section IV contains the results of the simulations and a comparison with the analytical predictions. II. ANALYTICAL THEORIES A. MSA for ions in a mixture of hard sphere solvents The description of the dielectric properties of solvents within the mean spherical approximation has been carried out by several authors starting from the work of Wertheim.23,24 Adelman and Deutch25 have solved the problem of a mixture of dipolar hard spheres in the framework of the MSA. In the case of equally sized components, they show that the mixture behaves like a monocomponent fluid with the same hard sphere radius and temperature parameter as the original mixture, but with effective dipole moment ␮ ˜ 2⫽1 m兺 k ␮ k 2,共1兲 and density ␳ ˜ ⫽1 ␮ ˜ 2兺 k ␮ k 2 ␳ k,共2兲 where mis the number of components in the mixture, and ␮ k and ␳ kare the dipole moments and number densities of each component, respectively. The correlation functions are completely determined by the correlation functions of Wertheim’s solution for the effective pure polar fluid. Chan et al.1have solved the problem of a mixture of ions modeled by hard spheres with point charges in a dipolar solvent of hard spheres with point dipoles, obtaining an expression for the free energy of solvation 共Born free energy of charging兲of an ion of charge ze at infinite dilution, FB⫽⫺ 共ze兲2 2共1 2R1⫹Rs兲 冉 1⫺1 ⑀ 冊 .共3兲 Here, ⑀ is the dielectric constant of the pure dipolar fluid, R1 is the ion diameter, and Rsis a correction to Born’s expression that is given in terms of the dipole radius and the dielectric properties of the pure dipolar fluid. Our aim is to determine the solvation energy of an ion in a binary mixture of dipolar fluids with equal radii but different dipolar moments. The result of Adelman and Deutch for the solvent mixture in the absence of ions suggests that Eq. 共3兲could be applied, using the dielectric properties of the effective fluid, as given by Eqs. 共1兲and 共2兲. In what follows, we will outline the necessary steps to generalize the theory of Chan et al. for a solvent consisting of two components with equal radii and different polarity. We will follow the notation of Chan et al. and will refer to their paper where appropriate. The ionic solutes are treated as charged hard spheres with diameter R1; different species of ions carry charges z ␣ and the electroneutrality condition 2361J. Chem. Phys., Vol. 113, No. 6, 8 August 2000 Solvation in dipolar binary mixtures 兺 ␣ ␳ ␣ z ␣ ⫽0共4兲 is assumed to hold, where ␳ ␣ denotes the number density for species ␣ . The ions are solvated by a mixture of hard spheres, each species having diameter R2, number density ␳ i, and carrying a dipole moment ␮ i. The interaction potentials have a hard core repulsive term and u ␣␤ 共r兲⫽z ␣ z ␤ e2/r,r⬎R1,共5兲 u ␣ i共r, ␻ 2兲⫽⫺z ␣ e ␮ iE2/r2,r⬎R12⫽1 2共R1⫹R2兲,共6兲 uij共 ␻ 1,r, ␻ 2兲⫽⫺ ␮ i ␮ jD12 /r3,r⬎R2,共7兲 where we have adopted Greek indices for ions and Latin indices for dipoles. The functions Ei⫽ ␮ ˆ( ␻ i)•rˆ and Dij ⫽ ␮ ˆ( ␻ i)(3rˆrˆ⫺I) ␮ ˆ( ␻ j) are the angular dependent parts of the interaction potentials. In these expressions, eis the elementary charge, Iis the 3⫻3 unit tensor, and the unit vector rˆ points from the molecule denoted by the first index towards the one denoted by the second index. The orientation of a dipole moment is defined by the solid angle ␻ . The Ornstein–Zernike 共OZ兲equations for an ion–dipole mixture are given by h ␣␤ 共r兲⫽c ␣␤ 共r兲⫹兺 ␥ ␳ ␥ 冕 dsc ␣␥ 共 兩 r⫺s 兩 兲h ␥␤ 共s兲 ⫹兺 k ␳ k 冕 ds 具 c ␣ k共r⫺s兲hk ␤ 共 ␻ 3,s兲 典 ␻ 3,共8兲 h ␣ j共r, ␻ 2兲⫽c ␣ j共r, ␻ 2兲⫹兺 ␥ ␳ ␥ 冕 dsc ␣␥ ⫻共 兩 r⫺s 兩 兲h ␥ j共s, ␻ 2兲⫹兺 k ␳ k 冕 ds ⫻ 具 c ␣ k共r⫺s, ␻ 3兲hkj共 ␻ 3,s, ␻ 2兲 典 ␻ 3,共9兲 hi ␤ 共 ␻ 1,r兲⫽ci ␤ 共 ␻ 1,r兲⫹兺 ␥ ␳ ␥ 冕 dsci ␥ ⫻共 ␻ 1,r⫺s兲h ␥␤ 共s兲⫹兺 k ␳ k 冕 ds ⫻ 具 cik共 ␻ 1,r⫺s, ␻ 3兲hk ␤ 共 ␻ 3,s兲 典 ␻ 3,共10兲 hij共 ␻ 1,r, ␻ 2兲⫽cij共 ␻ 1,r, ␻ 2兲 ⫹兺 ␥ ␳ ␥ 冕 dsci ␥ 共 ␻ 1,r⫺s兲h ␥ j共s, ␻ 2兲 ⫹兺 k ␳ k 冕 ds 具 cik共 ␻ i,r⫺s, ␻ 3兲hkj ⫻共 ␻ 3,s, ␻ 2兲 典 ␻ 3,共11兲 where his the total and cthe direct correlation function. We have used the notation 具典 ␻ ⫽1/4 ␲ 兰 d ␻ for the angular convolution. In the MSA, the closure conditions are obtained from Eqs. 共5兲–共7兲, hIJ共 ␻ 1,r, ␻ 2兲⫽⫺1, r⬍RIJ c, cIJ共 ␻ 1,r, ␻ 2兲⫽⫺ ␤ uIJ共 ␻ 1,r, ␻ 2兲r⬎RIJ c;共12兲 I⫽ ␣ ,i;J⫽ ␤ ,j. Here we have introduced upper case Latin indices which run over both Greek and Latin indices. In this notation the arguments of the correlation functions have to be adjusted depending on the current index: if two Greek indices are present, the argument is r. A Greek and a Latin index carry r and an angular dependence. In this case, if the Latin index is the second one, this angular dependence is ␻ 2, otherwise it is ␻ 1. Two Latin indices carry the full argument as in Eq. 共12兲. The contact distance RIJ cis R1for two Greek indices, R2for two Latin indices, and R12 for one Greek and one Latin index. The angular functions 1, Ei,Dij, and ⌬ij⫽ ␮ ˆ( ␻ i) • ␮ ˆ( ␻ j) form a closed set under the angular convolution 具 A( ␻ 1, ␻ 3)B( ␻ 3, ␻ 2) 典 ␻ 3. Multiplication tables for these quantities can be found in Ref. 1. Using this angular decomposition, we use an ansatz of the form f ␣␤ 共r兲⫽f11共r兲⫹z ␣ z ␤ fC共r兲, f ␣ j共r, ␻ 2兲⫽f12共r兲⫹z ␣ ␮ ˜ jfE共r兲E2, 共13兲 fi ␤ 共 ␻ i,r兲⫽f21共r兲⫺z ␤ ␮ ˜ ifE共r兲E1, fij共 ␻ i,r, ␻ 2兲⫽f22共r兲⫹ ␮ ˜ i ␮ ˜ jf⌬共r兲⌬12⫹ ␮ ˜ i ␮ ˜ jfD共r兲D12 for the correlation functions, where f⫽c,h, and ␮ ˜ i ⫽ ␮ i/ ␮ ˜ . This conveniently introduces the dipole moment of the effective fluid as defined in Eq. 共1兲. Note that our ansatz is essentially the same as the one in Ref. 1; the angular parts of the correlation functions have been scaled according to the different interaction potentials. The closure conditions 共12兲 in terms of these functions are the hard core conditions h11 ⫽h12⫽h21⫽h22⫽⫺1 for distances shorter than contact, and cC共r兲⫽⫺ ␤ e2 r,r⬎R1,共14兲 cE共r兲⫽ ␤␮ ˜ e2 r2,r⬎R12 ,共15兲 cD共r兲⫽ ␤␮ ˜ 21 r3,r⬎R2.共16兲 These conditions are the same as in Ref. 1, when ␮ in Ref. 1 is replaced by the dipole moment of the effective pure fluid ␮ ˜ . The three-dimensional Fourier transform f ˜ 共 ␻ 1,k, ␻ 2兲⫽ 冕 dreik•rf共 ␻ 1,r, ␻ 2兲共17兲 of the OZ equations in our notation may be written as h ˜ IJ共 ␻ 1,k, ␻ 2兲 ⫽c ˜ IJ共 ␻ 1,k, ␻ 2兲 ⫹兺 K ␳ K 具 c ˜ IK共 ␻ 1,k, ␻ 3兲h ˜ KJ共 ␻ 3,k, ␻ 2兲 典 ␻ 3.共18兲 The angular decomposition of f ˜ IJ( ␻ 1,k, ␻ 2) is the same as in Eq. 共13兲, when ris replaced by kand fby f ˜ . The 2362 J. Chem. Phys., Vol. 113, No. 6, 8 August 2000 Morillo et al. k-dependent angular functions will be denoted by 1, E ˜ i,⌬ ˜ ij, and D ˜ ij. Note that the coefficients of these angular functions may be obtained by one-dimensional Hankel transforms from their r-dependent counterparts f ˜ 共k兲⫽4 ␲ 共⫺i兲n 冕 0 ⬁dr jn共kr兲f共r兲,共19兲 where jn(x) is the spherical Bessel function of order nand the coefficients are transformed according to f ˜ 共k兲⫽f ˜ D共k兲;n⫽2, 共20兲 f ˜ 共k兲⫽f ˜ E共k兲;n⫽1, 共21兲 f ˜ 共k兲⫽f ˜ C共k兲,f ˜ ⌬共k兲,f ˜ ␣␤ 共 ␣ , ␤ ⫽1,2兲;n⫽0. 共22兲 In order to obtain decoupled equations for these quantities, it is convenient to introduce two new orthogonal angular functions J ˜ ⫹⫽1/3共⌬ ˜ ⫹D ˜ 兲,共23兲 J ˜ ⫺⫽1/3共2⌬ ˜ ⫺D ˜ 兲,共24兲 with coefficients f ˜ ⫹⫽f ˜ ⌬⫹2f ˜ D,共25兲 f ˜ ⫺⫽f ˜ ⌬⫺f ˜ D.共26兲 Carrying out the angular convolution in Eq. 共18兲, using the electroneutrality condition 共4兲, and setting ␳ c⫽兺 ␥ ␳ ␥ z ␥ 2,共27兲 ␳ I⫽兺 ␥ ␳ ␥ ,共28兲 ␳ d⫽兺 k ␳ k,共29兲 we obtain the following coupled set of equations for the correlation functions: h ˜ 11共k兲⫽c ˜ 11共k兲⫹ ␳ Ic ˜ 11共k兲h ˜ 11共k兲⫹ ␳ dc ˜ 12共k兲h ˜ 21共k兲,共30兲 h ˜ 12共k兲⫽c ˜ 12共k兲⫹ ␳ Ic ˜ 11共k兲h ˜ 12共k兲⫹ ␳ dc ˜ 12共k兲h ˜ 22共k兲,共31兲 h ˜ 21共k兲⫽c ˜ 21共k兲⫹ ␳ Ic ˜ 21共k兲h ˜ 11共k兲⫹ ␳ dc ˜ 22共k兲h ˜ 21共k兲,共32兲 h ˜ 22共k兲⫽c ˜ 22共k兲⫹ ␳ Ic ˜ 21共k兲h ˜ 12共k兲⫹ ␳ dc ˜ 22共k兲h ˜ 22共k兲,共33兲 h ˜ C共k兲⫽c ˜ C共k兲⫹ ␳ cc ˜ C共k兲h ˜ C共k兲⫺1 3 ␳ ˜ c ˜ E共k兲h ˜ E共k兲,共34兲 h ˜ E共k兲⫽c ˜ E共k兲⫹ ␳ cc ˜ C共k兲h ˜ E共k兲⫹1 3 ␳ ˜ c ˜ E共k兲h ˜ ⫹共k兲,共35兲 h ˜ ⫹共k兲⫽c ˜ ⫹共k兲⫺ ␳ cc ˜ E共k兲h ˜ E共k兲⫹1 3 ␳ ˜ c ˜ ⫹共k兲h ˜ ⫹共k兲,共36兲 h ˜ ⫺共k兲⫽c ˜ ⫺共k兲⫹1 3 ␳ ˜ c ˜ ⫺共k兲h ˜ ⫺共k兲.共37兲 This set of equations can be divided into three uncoupled groups: group A, Eqs. 共30兲–共33兲, together with the corresponding closure conditions 共12兲defines a mixture of hard spheres of diameter R1at density ␳ Iand diameter R2at density ␳ dwith the Percus–Yevick closure. Group B, Eqs. 共35兲–共36兲, describes the angular correlation of the fluid. Finally, Eq. 共37兲is decoupled from the other equations. When comparing our equations with the results of Chan et al.,we note that the angular correlations 共group B兲are described by the same OZ equations, when ␳ is replaced by the density of the effective fluid ␳ ˜ . As has been noted above, the closure conditions are also identical, replacing ␮ by ␮ ˜ . Thus, in what follows, we will use the results of Chan et al. for h ˜ ⫹ and h ˜ E. The average interaction energy of a single ion with the polar mixture is given by E⫽1 4 ␲ 兺 k ␳ k 冕冕 drd ␻ u ␣ k共r, ␻ 兲g ␣ k共r, ␻ 兲 ⫽⫺ 4 ␲ z ␣ 2e 3 ␳ ˜ ␮ ˜ 冕 R12 ⬁hE共r兲,共38兲 where we have used Eqs. 共6兲and 共13兲. The Born free energy of solvation for an ion of charge ze is then obtained as1 FB⫽⫺ 1 2 共ze兲2 1 2R1⫹Rs 冉 1⫺1 ⑀ 冊 , Rs⫽R2 冉 1 2⫺3 ␰ ˜ 共1⫹4 ␰ ˜ 兲 冊 ,共39兲 where ␰ ˜ is the solution of Q共2 ␰ ˜ 兲⫺Q共⫺ ␰ ˜ 兲⫽34 ␲ 9 ␤ • ␳ ˜ ␮ ˜ 2.共40兲 Here, Q( ␩ ) is the Percus–Yevick hard sphere inverse compressibility and ⑀ is the dielectric constant of the pure effective fluid25 共see the Appendix兲. In the following sections, we will compare the predictions of Eq. 共39兲with the simulation results. The solvation energy basically depends on the solvent longitudinal polarization density around the ion. The solvent polarization characterizes the response of the solvent to the field of the solvated ion. By definition, the polarization density of the kspecies in the mixture induced by the ion is given by Pk共r兲⫽ ␳ k 4 ␲ 冕 d ␻ g ␣ k共r, ␻ 兲 ␮ k共 ␻ 兲•rˆ⫽ ␳ k ␮ k 2z 3 ␮ ˜ hE共r兲, 共41兲 and the total polarization density is P共r兲⫽兺 kPk⫽1 3 ␳ ˜ ␮ ˜ zhE共r兲.共42兲 Far from the ion, the polarization density takes its asymptotic value given by the macroscopic expression Pmac⫽ze r2 1 4 ␲ 冉 1⫺1 ⑀ 冊 .共43兲 It will be illustrative to inspect the ratio Pk(r)/Pmac(r). This quantity represents the deviation of the polarization density with respect to the macroscopic model, and it empha- 2363J. Chem. Phys., Vol. 113, No. 6, 8 August 2000 Solvation in dipolar binary mixtures sizes the relative contribution of each species. The evaluation of the function hE(r) has been carried out in Ref. 1 and it leads to Pk共r兲 Pmac共r兲⫽ ␳ k ␮ k 2 ␳ ˜ ␮ ˜ 2 冉 1⫺rdF共r兲 dr ⫹F共r兲 冊 ,共44兲 where the quantity F(r) is the solution of an integral equation which can be generated numerically by a trapezoidal rule. The total polarization density can then be written as P共r兲 Pmac共r兲⫽1⫺rdF共r兲 dr ⫹F共r兲.共45兲 The relative contribution of species kis scaled by the factor ␳ k ␮ k 2. For the special case of a binary mixture with ␳ 1 ␮ 1 2/ ␳ 2 ␮ 2 2⫽1, both species contribute equally to the total polarization density P(r). These results will be used later on to discuss the meaning of preferential solvation of an ion in a mixture of fluids formed by hard sphere dipoles. B. The theory of Chandra and Bagchi Chandra and Bagchi2derived an expression for the energy of solvation of an ion in a binary mixture of polar liquids. Their work is based on an expansion of the free-energy functional of the inhomogeneous dipolar mixture in the presence of an external field 共the ion field兲. For the longitudinal polarization of the two species, they obtain Pi共k兲⫽3 4 ␲ ␥ i共k兲E0共k兲,i,j⫽1,2, 共46兲 where E0(k)⫽E0(k)k ˆis the bare electric field of the ion, E0(k)⫽i4 ␲ zej0(kR)/k, and ␥ i共k兲⫽YiiXjj共k兲⫺Xij共k兲Yjj Xii共k兲Xjj共k兲⫺Xij共k兲Xji ,i⫽j,共47兲 Xij共k兲⫽ ␦ ij⫺ ␮ i ␮ j ␳ i 4 ␲ cij m共110,k兲,共48兲 Yij⫽4 ␲ ␤ ␳ i ␮ i ␮ j/9. 共49兲 Here, the cij m(110,k) are coefficients in the expansion of the wave vector-dependent solvent–solvent direct correlation function in terms of spherical harmonics in the molecular frame 共this choice of frame is indicated by the superscript m兲.26 Namely, f共k, ␻ 1 ⬘, ␻ 2 ⬘兲⫽兺 l1l2mfm共l1l2m,k兲Yl1m共 ␻ 1 ⬘兲Yl2m ¯ 共 ␻ 2 ⬘兲, m ˜ ⫽⫺m.共50兲 The direct correlation function corresponds to the homogeneous dipolar mixture in the absence of solute ions and it may be obtained in analytical form from the MSA solution of the problem for the effective liquid.25 Thus, the coefficients cij m(110,k) can be readily evaluated within the MSA in terms of the Percus–Yevick solution for a system of nonpolar hard spheres with an effective density which depends on the polarity of the actual solvent. In what follows, we will restrict ourselves to a mixture of equally sized dipoles with diameter R⫽R2. For this case, the result of Chandra and Bagchi for the solvation energy is FB⫽⫺ 1 2 1 共2 ␲ 兲3 冕 dkP共k兲•E0共k兲 ⫽3共ze兲2 ␲ 冕 dk共 ␥ 1共k兲⫹ ␥ 2共k兲兲j0 2共kR兲,共51兲 and the polarization around the ion is given by P共r兲⫽P1共r兲⫹P2共r兲⫽6ze 共2 ␲ 兲2 冕 dkkj1共kr兲j0共kR兲 ⫻共 ␥ 1共k兲⫹ ␥ 2共k兲兲.共52兲 Using the definition of ␥ i(k) and Eq. 共A1兲, it is easy to show that P1共k兲 P2共k兲⫽ ␥ 1共k兲 ␥ 2共k兲⫽ ␳ 1 ␮ 1 2 ␳ 2 ␮ 2 2,共53兲 i.e., as in the previous approach, the ratio of the partial polarization densities is fixed by the system parameters. In order to numerically evaluate Eqs. 共51兲and 共52兲, the coefficients cij m(110,k) are required. In the Appendix, these functions are given in terms of the Percus–Yevick direct correlation function. The integrands in Eqs. 共51兲and 共52兲are difficult to treat numerically. It is convenient to study ␥ 12(k)⫽ ␥ 1(k)⫹ ␥ 2(k) for the limiting cases k→0 and k →⬁共see the Appendix兲and to carry out the integration in three parts, i.e., FB⫽3共ze兲2 ␲ 冋 ␥ 12共0兲 冕 0 k1dk j0 2共kR兲 ⫹ 冕 k1 k2dk ␥ 12共k兲j0 2共kR兲⫹ ␥ 12共⬁兲 冕 k2 ⬁dk j0 2共kR兲 册 , 共54兲 where the parameters k1and k2are chosen such that the result does not depend on these parameters. The integrals over j0 2(kR2) can be expressed in terms of sine integrals and the numerical integration can now be restricted to a finite range of values kwhere the integrands are numerically well behaved. III. MONTE CARLO SIMULATIONS We now turn to the calculation of the solvation free energy from Monte Carlo 共MC兲simulations of the model system considered in the previous section. Periodic boundary conditions with the minimum image convention27,28 were applied to a cubic simulation box of side length L. Convergence of the algorithm is accelerated by swapping solvent particles of different species. Details of our simulation method are described in Ref. 20. The simulations presented in this paper have been carried out using N⫽863 solvent particles. We used a packing fraction ␩ ⫽0.42, corresponding to a dense liquid. A cutoff of 4 molecular diameters has been applied for the long-ranged dipolar forces; contributions from outside the cutoff sphere have been taken into 2364 J. Chem. Phys., Vol. 113, No. 6, 8 August 2000 Morillo et al. account using a generalized reaction field method.29 For each system, 108MC configurations were generated; mean values were calculated after 107equilibration configurations. The dielectric constant of the reaction field has been obtained from simulations of the pure solvent. The free energy of charging of the solute ion in the solution is given formally by the thermodynamic integration expression30,31 ⌬F共0→1兲⫽ 冕 0 1d␭ 具 V 典 ␭,共55兲 which describes the free-energy difference between the uncharged state of the solute (␭⫽0) and its final charged state with total charge ze(␭⫽1). Vrepresents the total solute– solvent interaction energy and the angular brackets indicate an average taken with an equilibrium ensemble describing the overall system when the solute charge takes the fictitious value ␭ze. The numerical evaluation of the above expression requires, in principle, the knowledge of the integrand for a large number of points, thus requiring a substantial number of intermediate simulations. As pointed out by Hummer and Szabo,31 another alternative to evaluate the free-energy difference is to make use of the information about the statistical distribution of the values of Vavailable from computer simulations of just the initial and final states. Following their analysis, we evaluate the free-energy difference using the expression ⌬F共0→1兲⬇1 2共 具 V 典 0⫹ 具 V 典 1兲⫺ ␤ 12共 具 共V⫺ 具 V 典 0兲2 典 0 ⫺ 具 共V⫺ 具 V 典 1兲2 典 1), 共56兲 which depends on the knowledge of the first moments of the statistical distribution of Vin the initial and final states. As shown in Ref. 31, this expression is exact to fourth order in the free-energy perturbation expansion. In order to check the accuracy of Eq. 共56兲, we have carried out simulations for intermediate points; in all cases the corrections to the free energy are found to be within the statistical error limits. IV. RESULTS AND DISCUSSION In Fig. 1 we plot the Born free energy of solvation of an ion with valence z⫽1 as a function of the molar composition for a mixture of two hard sphere dipole liquids with equal radii and dipole moments ␮ L⫽0.17eÅ and ␮ H⫽0.34eÅ. The radius of the solute ion is taken to be the same as the solvent radius, R1/2⫽R2/2⫽1.44Å. Even for the cases of a single component solvent 共fH ⫽0 and fH⫽1兲, the solvation energy obtained from the MSA and the Chandra and Bagchi theory differs from the value obtained in the simulations. The deviation of the theoretical results from the simulation data is larger for the pure solvent with higher polarity (fH⫽1). The simulation results show two main features: first, a steep increase of FBcan be noted when a small fraction of high polarity solvent is added to the pure low polarity solvent (fH⬍0.1). Second, the value of the free energy of solvation becomes practically independent of the molar fraction for fH⬎0.5. This second feature of the solvation energy is better described by the MSA results and seems to be absent in the Chandra and Bagchi approach. Both theories reproduce the deviation of FB(fH) from an ‘‘ideal’’ law which implies a linear behavior with the molar fraction fH. This deviation, as measured by the excess free energy ⌬FB共fH兲⫽FB共fH兲⫺FB共0兲⫺fH共FB共1兲⫺FB共0兲兲,共57兲 has been used in Ref. 2 as an indication for preferential solvation. In Fig. 2 we plot ⌬FB(fH) for the same parameters as in Fig. 1. In this representation the deviation from the ideal behavior can clearly be seen for both the theoretical and the simulation results. However, the drastic increase of FBfor small molar fractions fH共which gives rise to a pronounced asymmetry in ⌬FB兲in the simulation results is not reproduced by the theories. The MSA approach yields a slightly higher value of the excess free energy when compared with the results of Chandra and Bagchi, although both theories underestimate ⌬FB. To further investigate the differences between the simulation results and the theoretical predictions, we now turn our attention to the induced polarization density P(r) around the FIG. 1. The Born free energy of solvation for the parameter values given in the main text. The circles are the simulation results; the solid line is an interpolation spline for these points. The dashed line is the theoretical prediction of Eq. 共39兲, the dotted-dashed line is the numerical result obtained from the Chandra and Bagchi approach, Eq. 共54兲. The long dashed line represents the macroscopic Born formula. The variable fHis the molar fraction of the component with higher polarity. FIG. 2. The excess free energy, Eq. 共57兲, for the same parameters and line types as in Fig. 1. 2365J. Chem. Phys., Vol. 113, No. 6, 8 August 2000 Solvation in dipolar binary mixtures ion. In Fig. 3 we represent the partial polarization densities PL(r)/Pmac(r) and PH(r)/Pmac(r) for the case of a mixture with ␳ L ␮ L 2/ ␳ H ␮ H 2⫽1, i.e., fH⫽0.2. For this choice of parameters, the partial polarizations of the low and high polarity species 关PL(r) and PH(r)兴obtained within the MSA 关see Eq. 共44兲兴 are identical. The same holds true for the Chandra and Bagchi approach 关see Eq. 共53兲兴. In Fig. 3, this is manifested by the coincidence of the two partial polarization density plots for each theoretical description. Both theoretical results give a similar behavior, except for values of rclose to contact. The first maximum of Pk(r)(k⫽L,H) lies at r ⫽Rin the MSA method, while the Chandra and Bagchi result shows a first maximum for values rslightly larger than R. For values r→⬁, both theories predict the expected macroscopic bulk value of the polarization density PL(r)/Pmac(r)⫽PH(r)/Pmac(r)⫽1/2. The failure of the theoretical polarizations to describe the simulation results are clearly seen in the figure. The partial polarization densities Pk sim(r) show oscillations, even for intermolecular distances at which the theoretical results have already reached their asymptotic value. Furthermore, the positions of the maxima and minima are displaced. Thus, the structure of the first few solvation shells is not predicted well by the theories. Another discrepancy manifests itself when studying the first and second solvation shell. While the theories predict PL(r)⫽PH(r) for all r, the simulation indicates that the contribution of PHprevails for small intermolecular distances. This feature is more pronounced for the first solvation shell. An interesting conclusion that can be drawn from the results discussed above concerns the meaning of preferential solvation. In the work of Chandra and Bagchi, preferential solvation is defined in terms of the deviation of the solvation energy from its ideal behavior with respect to the molar fraction. Let us suppose for a moment that the ion–dipole pair distribution function for species k⫽L,His defined as g ␣ k共r, ␻ 2兲⫽g ␣ k共r兲⫹z ␮ k ␮ ˜ hk E共r兲E2 k,共58兲 i.e., both the radial and the angular component of the correlation function depend on the index k. In this case, the angular average of ␮ k•rˆ is given by 具 ␮ kE2 k 典 共r兲⫽ 1 4 ␲ 兰 d ␻ 2 ␮ kE2 kg ␣ k共r, ␻ 2兲 1 4 ␲ 兰 d ␻ 2g ␣ k共r, ␻ 2兲 ⫽1 3z ␮ k 2 ␮ ˜ hk E共r兲g ␣ k共r兲,共59兲 and the partial polarization due to species kis Pk共r兲⫽ ␳ k 4 ␲ 冕 d ␻ 2 ␮ kE2 kg ␣ k共r, ␻ 2兲 ⫽ ␳ k ␮ k 2 具 E2 k/ ␮ k 典 共r兲g ␣ k共r兲⫽ ␳ k ␮ k 2h ˆk共r,fH兲,共60兲 where we have introduced h ˆk(r,fH). The argument fHin this function indicates that h ˆk(r,fH) generally depends on the molar composition of the mixture. The ratio of the partial polarizations will thus be given by PH共r兲 PL共r兲 ⫽ ␳ H ␮ H 2h ˆH共r,fH兲 ␳ L ␮ L 2h ˆL共r,fH兲 .共61兲 By writing the partial polarization densities in the form of Eq. 共60兲, we have conveniently separated the r-dependent contributions to Pk(r) that describe the angular correlations 具 E2 k/ ␮ k 典 (r), on one hand, and the radial distribution functions g ␣ k(r) on the other. In the previous sections, we have shown that for the two theories discussed in this paper, the ratio PH(r)/PL(r)isa constant independent of r:PH(r)/PL(r)⫽ ␳ h ␮ H 2/ ␳ L ␮ L 2. This arises because in the framework of the two theoretical treatments mentioned in this paper, h ˆ(r,fH) has the same functional form for both components of the mixture at a fixed molar fraction, so that h ˆH(r,fH)⫽h ˆL(r,fH)⫽h ˆ(r,fH). Note that the function h ˆ(r,fH) that is implicit in the work of Chandra and Bagchi differs from h ˆ(r,fH) in the MSA. The total polarization in terms of h ˆ(r,fH) is then P共r兲⫽PL共r兲⫹PH共r兲⫽ ␳ d共 ␮ L 2⫹fH共 ␮ H 2⫺ ␮ L 2兲兲f共r,fH兲. 共62兲 If h ˆ(r,fH) were independent of fH, Eq. 共62兲would yield a linear 共ideal兲behavior of the solvation energy with the molar fraction. Thus, the nonideality observed in Fig. 2 is a consequence of the dependence of h ˆ(r,fH) on the molar fraction. One should note that, even when using the Born formula for the solvation energy, with a dielectric constant of the mixture evaluated following the prescription of Adelman and Deutch,25 a nonideal behavior is obtained. This is to be expected, as the dielectric constant is not a linear function of fH. This limited description is not capable of describing the intuitive picture of preferential solvation: when the ion is preferentially solvated by the more polar species, one usually has in mind that the first solvation shell is composed, to a FIG. 3. The partial polarization densities Pk(R)/Pmac(R) around the ion for the same system parameters as in Fig. 1. The reduced radius Ris defined as R⫽r/R1⫽r/R2. The composition of the mixture is adjusted by the relation ␳ L ␮ L 2/ ␳ H ␮ H 2⫽1, i.e., fH⫽0.2. The simulation results PH(R) and PL(R) are represented by a solid line and a dotted line, respectively. For ␳ L ␮ L 2/ ␳ H ␮ H 2⫽1, both theories give PL(R)⫽PH(R). The MSA result is plotted as a dashed line, while a dotted-dashed line is used for the Chandra and Bagchi approach. 2366 J. Chem. Phys., Vol. 113, No. 6, 8 August 2000 Morillo et al. large extent, of molecules of the more polar species. This has as a consequence that g ␣ H(r) and g ␣ L(r) are different, and the ratio h ˆH(r,fH)/h ˆL(r,fH) is no longer unity and will generally depend on r. The simulation results for the partial polarizations shown in Fig. 3 show just this behavior. The ratio h ˆH(r,fH)/h ˆL(r,fH) differs from unity in the first two solvation shells; for distances close to contact we obtain values as high as h ˆH(r,0.2)/h ˆL(r,0.2)⬇25. It is instructive to analyze the ratios 具 E2 H/ ␮ H 典 ⫻(r)/ 具 E2 L/ ␮ L 典 (r) and g ␣ H(r)/g ␣ L(r) separately. In Fig. 4 we have plotted these ratios from the simulation data for a molar fraction of fH⫽0.2. It can clearly be seen that the main contribution to the deviation of h ˆH(r,fH)/h ˆL(r,fH) from unity stems from the radial contribution, i.e., from a spatial reorganization of the first two solvation shells. The ratio g ␣ H(r)/g ␣ L(r) reaches values of g ␣ H(r)/g ␣ L(r)⬇45 for distances close to contact. This corresponds to the intuitive picture of preferential solvation that we have given above. The angular contribution of the more polar species, 具 E2 H/ ␮ H 典 (r), when defined in this way, actually is smaller than its counterpart corresponding to the less polar species. However, the angular correlation 具 E2 H 典 (r)⫽ 具 cos( ␮ ˆH•rˆ) 典 is larger than 具 E2 L 典 (r), as expected. The preferential solvation parameter as defined by Marcus21,22 is ␦ fH⫽fH loc⫺fH,共63兲 where fH loc is the local molar fraction of species Hin a sphere of radius Rloc around the ion. Note that this parameter is zero both for the macroscopic 共Born兲and the MSA description as used in this work. In Fig. 5 we have represented this excess molar fraction for Rloc⫽4.0Å, corresponding approximately to the first minimum of g ␣ H(r). The simulation data show a quite large deviation from zero for a wide range of molar fractions. The maximum value of ␦ fH⬇0.6 is reached for fH⬇0.12, in agreement with the steep increase of the solvation free energy for small molar fractions, as has been noted above. The simulations indicate that a better theoretical description of the solvation of hard sphere ions in binary mixtures of hard sphere dipoles must take into account the rdependence of h ˆH(r,fH)/h ˆL(r,fH). Namely, an ion–dipole pair distribution function with the form of Eq. 共58兲, which contains separate radial and angular correlation functions for each species, should be used in order to adequately describe the phenomenon of preferential solvation. We will address this question in future work. The approach used here 关see Eq. 共13兲兴,as well as the theory of Chandra and Bagchi, do not fulfill this requirement. Nevertheless, they reproduce the qualitative behavior of the solvation energy with respect to the molar fraction and represent an improvement over the macroscopic description. ACKNOWLEDGMENTS Support by the Direccio ´n General de Ensen ˜ anza Superior of Spain 共Project Nos. PB98-0423 and PB98-1120兲and the Consejeria de Education y Ciencia of the Junta de Andalucı ´a is gratefully acknowledged. APPENDIX: MSA RESULTS FOR THE HOMOGENEOUS MIXTURE In the paper of Adelman and Deutch,25 the direct correlation function for a homogeneous dipole mixture of equally sized molecules of diameter ␴ is given in terms of the MSA solution for the effective fluid cij m共110,k兲⫽ ␮ ˜ ⫺2 ␮ i ␮ jcMSA m共110,k, ␳ ˜ , ␮ ˜ 兲,共A1兲 cMSA m共110,k兲⫽2共4 ␲ 兲2 ␬ 冕 0 ␴ drr2j0共kr兲cPY共r,2 ␬ ␳ ˜ 兲, 共A2兲 cPY共r兲⫽c0⫹c1 冉 r ␴ 冊 ⫹c3 冉 r ␴ 冊 3 .共A3兲 Here, cPY is the Perkus–Yevick direct correlation function at density 2 ␬ ␳ ˜ , and ␬ is obtained from the relation Q共2 ␬ ␩ ˜ 兲⫺Q共⫺ ␬ ␩ ˜ 兲⫽3y ˜ ⫽3共4 ␲ /9兲 ␤ ␳ ˜ ␮ ˜ 2,共A4兲 with Q共x兲⫽共1⫹2x兲2 共1⫺x兲4.共A5兲 FIG. 4. The ratios 具 E2 H/ ␮ H 典 (R)/ 具 E2 L/ ␮ L 典 (R) and g ␣ H(R)/g ␣ L(R) as obtained from the simulation. The parameters are the same as in Fig. 3. FIG. 5. The preferential solvation parameter as defined in Eq. 共63兲. The local molar fraction has been determined from the radial distribution functions inside a sphere with radius Rloc⫽4.0 Å around the solute. The circles represent the simulation results; the solid line is a spline interpolation of these points. 2367J. Chem. Phys., Vol. 113, No. 6, 8 August 2000 Solvation in dipolar binary mixtures