Improvements of the continuum model. 1. Application to the calculation

J . Phys. Chem. 1988, 92, 1617-1631

1617

Improvements of the Continuum Model. 1. Application to the Calcuiatlon of the Vaporization Thermodynamic Quantities of Nonassociated Liquids J. Langlet,*+P. Claverie: J. Caillet,? and A. Pullman Institut de Biologie Physico-Chimique, Laboratoire de Biochimie Thgorique Associe au CNRS 13, Rue Pierre et Marie Curie, 75005 Paris, France, and Laboratoire de Dynamique des Interactions MolCculaires du CNRS, E R No. 271, Universite Pierre et Marie Curie, 4, Place Jussieu, 75230 Paris Cedex 05, France (Received: December 9, 1986; In Final Form: May 14, 1987)

A significant effort has teen made in order to improve the continuum model used for calculation of the solvation thermodynamic quantities of a molecule embedded in a cavity formed by the intersecting van der Waals spheres of the solute in a polarizable medium. These improvements principally concern the thermodynamic quantities associated with the electrostatic part of the solvation energy. A simple method is proposed for the calculation of the fictive charge density representing the reaction potential. The calculated values obtained agree quite well with those calculated within more sophisticated methods. The improvements also principally concern the representation of solvation sites. The method is applied to the calculation of the vaporization thermodynamic quantities of nonassociated liquids, and the results obtained are discussed in relation with experimental data.

I. Introduction In recent years, there has been a growing interest in the possibilities of theoretical interpretation of the effect of solvents on the properties of molecules (in particular biomolecules) in relation to conformations (ref 1-4 and references cited therein) and to complexing a b i l i t i e ~ . ~ - ~ ’ The theoretical study of solutions is made difficult by the very large number of molecules involved, with the difficulty proceeding from the problem of the simultaneous evaluation of their mutual interactions. Monte Carlo methods of statistical mechanics have been increasingly used recently,22-28providing information about ensemble averages of mechanical properties and solvent microstructures. Previously, Monte Carlo methods only dealt with the ~ - ~ ~effort has been solvation enthalpy PH, but r e ~ e n t l y a~ great made by using Monte Carlo methods and the coupling parameter approach to determine the solvation free energy AGsalvwhich is very important for comparing theoretical and experimental data since the equilibrium constants measured experimentally are related to AG. Such calculations, however, require considerable computing time, so there is still a need for simpler methodologies. In this area, aside from or rather complementary to what can be called “small clusters” method^,^^-^' it appears worth developing accurate “macroscopic” or “continuous” methods. The well-known principle of such rests on the representation of the solute molecule as a charge distribution placed inside a cavity surrounded by a continuous medium standing for the solvent characterized by macroscopic properties: dielectric constant, molecular volume, thermal expansion coefficient, etc. The aim is the determination of both free enthalpy AG and enthalpy AH (hence the entropy). Although such a model is very appealing due to its apparent simplicity, its practical implementation presents specific problems, in particular (a) the choice of the shape of the cavity containing the solute molecule, (b) the description of the solute-solvent interaction a t short distances, (c) the correct introduction of the electrostatic interactions, (d) the representation of the reorganization of the solvent in the presence of the solute molecule, and (e) the proper representation of the “solvation sites”, i.e., the parts of the solute molecule that interact very strongly with the solvent. Various solutions for each of these problems have been proposed in recent years (see for instance ref 39-41 and references cited therein), making possible the development of relatively accurate continuum models. One such model has been elaborated in our laboratory39on the basis of a methodology initially proposed for pure liquid hydrocarbons!2 We have now improved this model particularly in the part concerning the solutesolvent electrostatic interactions, as well as in the practical. Universite Pierre et Marie Curie.

computational procedures of the different terms, aiming a t a reasonable compromise between accuracy and applicability to large (1) Pullman, A.; Pullman, B. Q.Reu. Biophys. 1975, 7, 505. (2) Avignon, M.; Garrigou-Lagrange, C.; Bothorel, P. Biopolymers 1973, 12, 1661. (3) Ison, R. J. Pharmacol. 1972, 24, 82. Ison, R.; Partington, P.; Roberts, C. C. Mol. Pharmacol. 1973, 9, 958. (4) Kuntz, I. D.; Kauzmann, W. Adu. Protein Chem. 1974, 28, 239.

( 5 ) Kuecher, E.; Derkosch, J. Z . Naturforsch., B: Anorg. Chem., Org. Chem. 1966, 218, 209. (6) Hamlin, R.; Lord, R. C.; Rich, A. Science (Washington, D.C.) 1965, 148, 1734. (7) Pitha, J.; Jones, R. N.; Pithova, P. Can. J . Chem. 1966, 44, 1045. (8) Katz, L.; Penman, S. J. J . Mol. B i d . 1966, 15, 220. (9) Kyogoku, Y.; Lord, R. C.; Rich, A. J . Am. Chem. SOC.1967,89,496. Kyogoku, Y.; Lord, R. C.; Rich, A. Biochim. Biophys. Acta 1969, 179, 10. (10) Miller, J. H.; Sobell, H. M. J. J . Mol. Biol. 1967, 23, 345. (11) Binford, J. S.;Holloway, D. M. J . Mol. Biol. 1968, 31, 91. (12) Ts’o, P. 0. P.; Melvin, I. S.; Olson, A. C. J . Am. Chem. SOC.1963, 85, 1149. (13) Schweizer, M. P.; Broom, A. D.; Ts’o, P. 0. P.; Hollis, D. P. J . Am. Chem. SOC.1968, 90, 1042. (14) Ts’o, P. 0. P.; Kondo, N. S.; Robins, R. K.; Broom, A. D. J. Am. Chem. SOC.1969, 91, 5625. (15) Helmkamp, G. K.; Kondo, N. S.Biochim. Biophys. Acta 1967, 143, 27; Biochim. Biophys. Acta 1968, 157, 242. (16) Nakano, N. I.; Igarshi, S. J. Biochemistry 1970, 9, 577. (17) Evans, F. E.; Sarma, R. H. Biopolymers 1974, 13, 2117. (18) Imdo, T. Biochim. Biophys. Acra 1977, 475, 409. (19) Plisiewicz, E.; Stepien, E.; Wierzchowski, K. L. Nucleic Acids Res. 1976, 3, 1295. (20) Mayevsky, A. A.; Sukhorukov, B. I. Nucleic Acids Res. 1980,8,3029. (21) Farquhar, E. L.; Downing, M.; Gill, S.I. J . Biochem. 1968, 7, 1224. (22) Clementi, E.; Corongiu, G. J . Chem. Phys. 1980, 72, 3979. (23) Jorgensen, W. L. J. Am. Chem. Soc. 1979,101,2016. Jorgensen, W. L.; Binning, R. C.; Bigot, B. J. Am. Chem. SOC.1981,103,4393. Jorgensen, W. L. J. Am. Chem. SOC.1981,103,4721. Jorgensen, W. L. J. Chem. Phys. 1982, 77, 5757. Jorgensen, W. L. J. Am. Chem. SOC.1981, 103, 335, 341, 345. Jorgensen, W. L.; Ibrahim, M. J. J. Am. Chem. SOC.1981,103,3976. Chandrasekhar, J.; Jorgensen, W. L. J. Chem. Phys. 1983, 77, 5073. Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J . Am. Chem. SOC.1984,106, 6638. (24) Cloessens, M.; Ferrario, M.; Ryckaert, J. P. Mol. Phys. 1983,50, 217. (25) Linse, P. J. Am. Chem. SOC.1984, 106, 5425. (26) Beveridge, D. L.; Mezei, M.; Mehrotra, P. K.; Marchese, C. T.; Ravishanker, G.; Vesu, T. R.; Swaminathan, S. In Molecular Based Study

and Prediction offluidproperties; Haile, J. M.; Mansoori, G., Eds.; Advances in Chemistry 24; American Chemical Society: Washington, DC, 1983. Beveridge, D. L.; Mezei, M.: Mehrotra, P. K.; Marchese, F. T.;Thirumalai, V.; Ravishenker, G. Ann. N.Y. Acad. Sci. 1981, 367, 108. Mezei, M.; Mehrotra, P. K.; Beveridge, D. L. J . Biomol. Struct. Dyn. 1984, 2, 1. Beveridge, D. L.; Meye, P. V.; Jayrram, B.; Ravishanker, G.; Mezei, M. J . Biomol. Struct. Dyn. 1984, 2, 261. (27) Danilov, V. I.; Tolokh, I. S. J . Biomol. Struct. Dyn. 1984, I , 119. Danilov, V. I.; Tolokh, I. S.; Poltev, V. I.; Malenkov, G. G. FEES Lett. 1984, 167, 245. (28) Pohorille, A,; Pratt, F. R.; Burt, S. K.; Mac Elroy, R. D. J. Biomol. Struct. Dyn. 1984, I , 1257. Pohorille, A,; Burt, S. K.; Mac Elroy, R. D. J . Am. Chem. SOC.1984, 106, 402.

0022-3654/88/2092-1617$01,50/00 1988 American Chemical Society

Langlet et al.

1618 The Journal of Physical Chemistry, Vol. 92, No. 6, 1988

systems of biological interest. Essential to the present developments are (i) a reexamination of the derivation of the thermodynamical quantities associated with the electrostatic part of the solute-solvent interactions, (ii) the choice of the shape of the cavity as the van der Waals envelope of the solute molecule, and (iii) the representation of the solute within the cavity by a set of multipoles (up to the quadrupoles) derived from a b initio wave function. Furthermore, a procedure to take into account the existence of “solvation sites” is implemented by defining in a systematic way appropriately modified van der Waals radii for each atom of the solute molecule. In this paper, we present the resulting methodology and its application to the calculation of the vaporization energies (Le., the opposite of the solvation energy of vapor condensing into its liquid) of various pure liquids made of nonpolar and nonassociated polar molecules.

In eq l a A is the area of the boundary surface of the cavity. Consistent with our choice of the cavity shape we adopt for A the sum of the areas Ai of the fractions of spheres implicated. In eq 1b aD is the thermal expansion coefficient of the solvent and T the temperature in absolute degrees. (b) Since the surface created is of microscopic dimensions and is highly curved, it may be preferable43to use rather AEcav = xg(b/a)yA[ 1 - T($

+

$)]

(2a)

where K:(b/a) transforms the macroscopic surface tension y to microscopic dimensions by

11. Methodology

Let us recall that, in a continuum model, the calculation of the solvation energy involves (a) the calculation of the cavitation energy (energy necessary to create in the solvent the appropriate cavity of the solute molecule), and (b) the calculation of the energies corresponding to the introduction of the solute molecule into the cavity, namely, the thermodynamic quantities AG, AS, and AH (free enthalpy, entropy, and enthalpy) corresponding to the solvation process. (A) Cavitation Energy. The calculation of the cavitation energy necessitates the choice of the shape of the cavity: we adopt the ~ imchoice proposed by Huron and C l a ~ e r i e(a~ considerable provement over earlier representations), namely, the utilization of a volume corresponding to the exact molecular shape of the solute, in practice made of the union of the van der Waals spheres of the atoms composing it. Thus the cavitation energy is calculated as a sum of contributions from the corresponding spherical atomic fragments. For the practical calculation of these terms we shall compare the results of three techniques: the first one based on the classical derivation using the macroscopic surface tension y, the second one using a factor that transforms the macroscopic surface tension to microscopic dimension^:^ and the last one utilizing the Reiss-Pierrotti expressions(based on the “scaled particule” theory of for the energy of cavitation of spherical cavities. The respective expressions utilized are: (a) In terms of the surface tension y the free energy of cavitation AFav and the corresponding internal energy AE,, are given respectively by AFav = yA

(la)

~~

(29) Beveridge, D.L.;Mezei, M.; Ravishanker, G.; Jayarem, B. Int. J . Quantum Chem. 1986, 29, 1513. (30) Shing, K.S.; Gobbins, K. E. Mol. Phys. 1981, 43, 717; Mol. Phys. 1982, 46, 1109. (31) Mehrotra, P. K.;Mezei, M.; Beveridge, D. L. J . Chem. Phys. 1983, 78, 3156. (32) Mezei, M.Mol. Phys. 1982, 47, 1307. (33) Pullman, A. In Quantum theory of chemical reactions; Daudel, R.; Pullman, A.; Salem, L.; Veillard, A,, Eds.; Reidel: Dordrecht, 1980 Vol. 11,

v b and Va are the molecular volumes of the solvent and solute, respectively. When dealing with pure liquids, eq 2b only contains the term Kb( 1) which has been estimated from solubility Within this approach, the free energy of cavitation is calculated as

with

and

where K”,b/a) is a constant similar to &b/a) calculated by using xi( 1). This term occurs in the entropy term and measures the deviation of the cavity entropy from the macroscopic surface entropy.43 In both cases the free enthalpy of cavitation AGaV is obtained as AGaV = AFcav + PVa

(c) In the third p r o c e d ~ r e ,the ~ ~ cavitation *~~ free enthalpy AGF(R,) of a sphere of diameter 2Ri in a solvent considered as made itself of spheres of diameters 2a is AGF(Ri) = KO + K , d

+ K2dz + K3d3

(3)

where d = 2(Ri + a) is the diameter of a sphere that excludes the centers of the solvent molecules, and

P 1.

(34) Pullman, B.;Miertus, S.; Perahia, D. Theor. Chim. A r f a 1979, 90, 317. (35) Goldblum, A.; Perahia, D.; Pullman, A. FEBS Lett. 1978, 91, 213. (36) Mantione, M. J.; Daudey, J. P. Chem. Phys. Lett. 1970, 6 , 93. (37) Huron, M. J.; Claverie, P. Chem. Phys. Lett. 1971, 9 , 194. (38) Egan, J. T.;Nir, S.; Rein, R.; MacElroy, R. Int. J . Quantum Chem., Quantum Biol. Symp. 1978, 5, 433. (39) Claverie, P.; Daudey, J. P.; Langlet, J.; Pullman, B.; Piazzola, D.; Huron, M. J. J . Phys. Chem. 1978, 82, 405. (40) Miertus, M.; Tomasi, J. Chem. Phys. 1982,65, 239. (41) Claverie, P. Quantum Theory of Chemical Reactions; Daudel, R.; Pullman, A.; Salem, L.; Veillard. A,, Reidel: Dordrecht, 1980; Vol. 11, p 151. (42) Huron, M. J.; Claverie, P. J . Phys. Chem. 1972, 76, 2123. (43) Halicioglu, T.; Sinanoglu, 0. Ann. N.Y. Acad. Sei. 1969, 158, 308. (44) Reiss, M.; Frisch, H. L.; Helfand, E.; Lebowitz, J. L. J . Chem. Phys. 1960, 32, 119. (45) Pierotti, R. A. J . Phys. Chem. 1963, 67, 1840. (46) Pierotti, R. A. J . Phys. Chem. 1965, 69, 281.

(2e)

K2 = N E

4a2

[ 12y + 1-Y

18(

e)’]

K3 = K r P

- 2rPa

(4)

withy = 8 r a 3 ( p / 6 ) , p being the density number of the solvent, N Avogadro’s number, k the Boltzmann constant, P the pressure, and T the absolute temperature. Within the model adopted for the cavity, each fragmental sphere is taken as contributing a fraction of AG ( Ri ) proportional to its area Ai, so that (see ref 39)

The Journal of Physical Chemistry, Vol. 92, No. 6, 1988 1619

Improvements of the Continuum Model and the corresponding enthalpy can be deduced by the Helmholtz expression:

interaction can be written as (assuming for simplicity monopoles gj only on the solute) (C

According to Pierotti:5*4 eq 3 has the form required in classical thermodynamics. In effect, when neglecting the very small term involving P which is just the volume work, we notice in eq 2c a term involving the surface tension (depending on 62,which is a surface work) and a term depending on d which corrects the surface tension for the effect of surface curvature. As noticed by Pierotti, the only difference between eq 3 and those given by classical thermodynamics (eq 2c, for instance) is the absence of the constant term KOin the thermodynamic equation: this term is important for microscopic cavities. The results of the three above formulations will be compared in the present paper. We shall see in section V1.A that it is far from indifferent to use one or the other. ( B ) Solute-Solvent Energies. The energy changes accompanying the insertion of the solute into the cavity are best introduced by using the charging parameter method4' between two thermodynamic states 0 and 1 characterized by intermolecular energy functions U,and VI,respectively, the transformation from state 0 to state 1 being made with the help of a charging parameter f (varying from 0 to 1) and an energy function U([)such that

U(0) =

u,

U(1) =

u,

Hence, the solute-solvent interaction is suppressed when f = 0 and may be continuously "switched on" up to its real physical value by varying f from 0 to 1. The amount of work associated with such a reversible charging process corresponds to the variation AG of the free enthalpy between the uncharged and charged states of the system by the integral

= CqjJ.I(Fj)

s v )1

(8)

j

in terms of the reaction potential + I of the solvent fully polarized by the solute fully charged. The reaction potential is the statistical average of the electrostatic potential created by all the solvent molecules. Since the distribution function of the solvent molecules depends on the degree of charging of the solute, so does the reaction potential. Let us denote $&Fj) the reaction potential corresponding to the level [ of charging. The usual theory of ) respect to dielectrics assumes a linear dependence of I++(?,with

5: =

+E(?,)

(9)

t$l(Fj)

The free enthalpy variation Affl upon charging can be obtained as the reversible work of charging the solute. In the elementary step where the charging parameter increases from f to f d[, the charges vary from fq, to (f df)q,, hence dqj = qj df. These additional charges are brought into the reaction potential through the elementary work dw:

+

+

dw = C+E(Fj)qj df

(10)

dw = ?t$l(Fj)qj d t

(1 1)

J

or, according to (9)

where f$q(Yj)q, corresponds to the average ( dU/df)' appearing in eq I. Integrating over f from 0 to 1 yields

(12a) AGC1= YZFqj+lVj) over the canonical average of the derivative of the interaction energy U. In continuum models this canonical averaging is replaced by a macroscopic-type calculation which implies an appropriate choice of the initial and final states in such a way that the associated switched-on interactions correspond at best to an average interaction potential easy to calculate (see ref 39 and 41 for a discussion). In the process adopted here, the initial step is the creation of the cavity as the previously described van der Waals volume of the solute inside the solvent with bulk properties including all solvent interactions. Thus for the "charging" process, the initial state (f = 0) corresponds to the empty cavity surrounded by the solvent, the interactions of which are completely switched on. Then, when [ varies from 0 to 1, the solutesolvent interactions are progressively switched on, and in the final state (f = l), all the interactions between the solute and on the solvent are switched on. In this charging process the parameter f switches on only the solutesolvent interactions. For convenience, the charging process is subdivided into two parts involving respectively (a) the switching on of the electrostatic part and (b) the switching on of the dispersion-repulsion part of the interactions. We examine first, in the next section, the derivation of the expression of the free enthalpy involved in part a and its consequences and then the practical evaluation of the corresponding quantities within the framework of our representation. ( a ) The Electrostatic Part of the Solute-Solvent Interactions. According to the preceding remarks, the initial state of this part corresponds to the uncharged solute inside the cavity and the final state to the fully charged solute and to the solvent having undergone a reorganization under the polarizing effect of the solute. In this final state ( f = l), the solute interacts with the reorganized polarized solvent. The electrostatic part of the corresponding (47) Kirkwood, J. G. Theory ojLiquids; Alder, B. J., Ed.;Gordon and Breach: New York, 1968.

( 12b)

This expression of AGcl compared to eq 8 leads to the relation AG~I

=

y2(q!sv )

(13)

which indicates that the charging work is equal to half the solutesolvent corresponding electrostatic interaction energy. Note that the appearance of the coefficient 1/2 in relation 13 comes from the linear dependence of the reaction potential on f , a feature that was not quite clearly understood earlier (see, e.g., ref 39). In fact, contrary to the statement in ref 39, the AGC1calculated is the complete AGC1corresponding to switching on the solute-solvent electrostatic interaction. Equation 12b is readily extendable from charges to multiples of higher degrees, the corresponding parts of the reaction potential being then replaced by its derivatives which also show a linear dependence on E. In our model the charge distribution of the solute is represented by a b initio (within an appropriate basis set48) multipoles (up to quadrupoles) located on atoms and bonds, calculated by a procedure defined by VignbMaeder and Claverie." Coming back to relation 13 between AGel and ( C&v)l, let us show that it provides a way to obtain the electrostatic part of the solvent reorganization energy as follows: on the one hand, we can from Ace' obtain the corresponding enthalpy MI by the standard relation =

-T

a(AG"') aT

On the other hand, the internal energy variation is the difference and (VI) between ( I which are respectively: (48) Berthod, H.; Pullman, A. J . Comput. Chem. 1981, 2, 87 and references cited therein. (49) Vi@-Maeder, F.; Claverie, P. J . Chem. Phys., in press. Vi@Maeder, F.; Claverie, P. J . Mol. Struct., TheoChem. 1984, 107, 221.

1620 The Journal of Physical Chemistry, Vol. 92, No. 6, 1988 (%v)o

and

(esv),

+ (K;,,),

Langlet et al. u on the surface of the c a ~ i t y . ~ ~ ~ ~ ~

thus AW' = (G!sv)~ + AuSv,, (15) where AU,, is the solvent reorganization energy due to the switching on of the electrostatic part of the energy. This quantity can thus be obtained as

Aqi-sv

AHC'

- 2AGC'

(16)

using eq 15 and 13. In the above section, we have used to thermodynamical quantities AG" (hence AHC'). In fact, the procedure allows us to calculate either A@ or A P 1 (hence according to whether the elementary work dw is calculated at constant pressure P or at constant volume V. The justification of the use of AGC1in our methodology will be discussed at the end of this section. The practical calculation of AGel using eq 8-13 requires the evaluation of the reaction potential which is a statistical average of the electrostatic potential created by the solvent molecules. In practice, in a continuum model, J / , is evaluated approximately by the reaction electrostatic potential V' evaluated within the framework of a dielectric model (the solvent represented by a dielectric medium). We describe below how y e proceed. Under the influence of the electrostatic field Qo created by the charge distribution (multipolar one in our model) of the solute, the outside dielectric medium (the solvent) is polarized and acquires a dipole moment 8 per unit volum_e, which produces a reaction potential V'(thus a reaction field Q'), hqnce a resulting total potential V and a total electrostatic field 6. V=vO+V' (17) 8 = 20 + 2' (18)

e),

The dipole moment density 8 acquired by the solvent is given (using rationalized system of units) by d = (1/4lr)(t - to)$ (19) where to is the vacuum dielectric permittivity, t is the permittivity of the solvent, and €/eo is the usual dielectric constant tr. The potential Vcan be obtained by solving the Laplace equation AV=O (20) with the necessary boundary conditions on the cavity surface and a t infinity (see, for instance, ref 50 and 51): v(qlI = V(q)IE (21)

aV(q

aV(q

XI,= XIc

with

(22)

=0

(23) =0 (24) with 1, and IE indicating points lying at infinitely small distances from the surface, respectively inside and outside the cavity. The solution of the problem which is known in very simple cases, such as spherical or prolate spheroidal c a ~ i t y , becomes ~ ~ - ~ ~more difficult for a cavity of arbitrary shape. Huron and C l a ~ e r i e ~ ~ proposed a numerical method where the potential is developed in a basis of harmonic functions. This, however, requiring relatively lengthy computations, we turned to another procedure using the fact that, in an homogeneous dielectric medium (the dielectric constant e, = t/tOindependent of 7) which does not contain any real charge distribution, the reaction potential V'created by the dipolar density 9 is equivalent to that of a fictive charge density

Ilrl*lgrad

(50) Rocard, Y. ElectricirP; Masson: Paris, 1956. (51) Durand, E. Elecrostatique et MagnZtostatique; Masson: Paris, 1953. (52) Kirkwood, J. G. J . Chem. Phys. 1984, 2, 351. Kirkwood, J. G.; Weistheimer, F. M. J. Chem. Phys. 1936, 6, 506. (53) Onsager, L. J . Am. Chem. SOC.1938, 58, 1486. (54) Born, M. J. Phys. 1920, 1 , 45. (55) Weistheimer, F. M.; Kirkwood, J. G. J . Chem. Phys. 1938, 65, 13. (56) Abbott, J. A.; Bolton, H. C. Trans. Faraday SOC.1952, 48, 422. (57) Wada, A. J . Chem. Phys. 1954, 22, 198. (58) Buckingham, D. Trans. Faraday SOC.1953, 49, 881. (59) Huron, M.J.; Claverie, P. J . Phys. Chem. 1974, 78, 1853.

u

=

-

1 41r

- €o)z.Q

-(t

ii being the unit vector normal to the surface.

Equation 25 requires the calculation of integral equation5' € +

1

-u(?)

t-1

= 2to(Z,ho)

u which is given by

1 + -l$u(B') 21r

the

cos 0 7 ds" (28)

1s - S I *

where s' and s" denote two points on the cavity surface. The derivation of eq 28 is recalled in Appendix I, where we also briefly discuss the connection between this equation and the procedure of Miertus et a1.60 The calculation of u, using eq 28, requires in principle an iterative process. But this again requires relatively lengthy computations. In order to avoid this difficulty, we have used a simpler procedure taking advantage of our multipole representation of the charge distribution of the solute and of the representation of the cavity by the union of van der Waals atomic spheres. First of all, it must be emphasized t i a t , due to the linear the total density character of eq 28 with respect to u and Q04, u may be obtained as the sum of partial densities u(i,l), where u(i,l) denotes the density that would be generated by the multipole of order 1 located at the center i supposed to be the only one present in the cavity. Writing

ii-80= Ci C/z & ( i , l )

(29)

and inserting eq 29 in eq 28 and using the linear character of the latter, we actually obtain u

=

CCC(i,I) i /

(30)

where u(i,l) is the charge density created by i i . $ O ( i , l ) alone. The problem now reduces to express each partial density .(ill) in some simplified way (avoiding the exact iterative process). For that purpose, y e shalluse the facithat, in the special case of a spherical cavity, &'(I), Q'(l), hence e(/), and u(l) are exactly known (see Borns4for 1 = 0 (charge), Onsagers3for I = 1 (dipole), and Kirkw@** for a-multipole of any order I). Once the expressions of Qo(l)and Q'(l) are known, it is easy to calculate their Gormal compnent to the surface and to find a relation between go,,(l)and 6,,(Z), namely

&(I)

= f(l):on(l)

(31)

The factorf(1) has been calculated for a charge and a dipole.60 We have generalized the calculation to a multipole of any degree 1 (see Appendix II), yielding the result (21 + I)(€ - €0) = (I 1)t leo

+

+

Thus in terms of the normal component of the parts of 8 O due to the multipoles of different degree I , the charge surface density for a sphere embedded in a dielectric medium can be calculated exactly as (21

+ l ) ( e - eo)

t[ ( I + 1)e +

u = 1-0

It0

8°.(l)]

(33)

In Appendix I1 we show that the factor f(I) calculated by eq 32 is different from the factor given by Miertus et aL60 but that u remains the same in the two cases. (60) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117.


Improvements of the Continuum Model

Cedi =

f"

s,.(s)P(s)

ds'

(39)

Thus, the solute-solvent electrostatic energy appears as the energy of the fictive surface density u in the electrostatic potential is done by a surface created by P O ( ? ) . The calculation of (&,) integration using a Korobov grid as in ref 42, using for P(8the expressions derived from the utilization of our multipolar distribution: N

P(i)=

i- 1

+

P(i,O) ?"'(i,l)+ P ( i , 2 )

(40)

with

Figure 1. Geometric buildup of the component, perpendicular to the solute cavity surface at point P, octhe electric field due to the monopoles ( I = 0) of the solute molecule: &",(P) = ~ i ~ , i i . & o ( i , l , P ) .

The utilization of eq 3 3 in the case of a cavity formed by the union of van der Waals spheres (Figure 1 ) is done as follows: In our model with multipolar distribution centered on atoms and middle of bonds, each atom ( i = 1 , 3 in Figure 1) and each middle of a bond ( i = 4 , 5 in Figure 1) carries multipoles up to quadrupoles (I = 2 ) . We then calculate at any poi$ P of the cavity surface the normal components eo,(i,l) = X o ( i , 1 )of the electrostatic field created by the multipole of order I belonging to the center i, and we express in an approximate y a y the corresponding partial charge density u(i,l) in terms of eo,(i,l) by using eq 33 (which would be exact for a spherical cavity only), namely u(i,@) = f ( r ) &,(i,@)

(34)

hence the approximation for the total charge density at F N

a(?) =

C v ( 4 Zon(i,/;F) i=l I

(35)

Note that, even when 3 belongs to the sphere with center i, expression 34 cannot be strictly exact (beyond the special case of a cavity reduced to this single sphere), because the exact u(i,l) depends upon the shape of the full cavity, which is involved in the boundary conditions of the basic system of eq 21-24; we can only expect thaieq 34 gives a reasonably accurate approximation. The part of u(P) due to multipoles located at centers different from i may of course be reproduced still less accurately, but since these contributions are also smaller, because these centers are more distinct from P , the overall accuracy of the calculated total u can remain satisfactory. Once u has been calculated by eq 3 5 , we obtain V'using eq 25, and ( q L v ) , (eq 12b) may be written as (&)I

= sV'(9

PO(?)

d7

(36)

with PO(?) denoting the charge distribution of the solute (which in eq 12 was reduced to discrete charges q only). Equation 36 may be rewritten as

(43) where N is the number of charged centers (on atoms and bonds) and q, and 6, represent respectively the monopole and dipole on center i. The total quadrupole is reduced to two axial quadrupoles of magnitude 8, with direction a', (see ref 4 9 ) . The AGel and AF1are calculated according to eq 12b and 14, respectively. The derivative d/dT(AGC!) entering the evaluation of AF! is calculated in a similar way as in ref 59: AGCI depends on temperature T through the dielectric constant and the surface S limiting the cavity. Note that in our methodology, the calculation of AGe' using eq 12, 13, 32, 3 3 , and 39 and of AF1using eq 14 involves macroscopic parameters given at constant pressure P the dielectric constant is given at P = 1 atm and the derivatives (atr/aT), and (aV/dT), (the expansion coefficient) are given at P constant. Thus during the switching on of the electrostatic part of energy, the only neglected contribution dw of the work evaluated under constant pressure is PAV, where AV is the variation of volume of the solution due to the switching on of this interaction. Since AVis very weak, and the corresponding work is very weak kcal/mol for (Ben Naim61v62has evaluated PAV N 1.2 X the dissolution of methane in water), we are justified in neglecting this contribution. (b) Expression of the Solute Polarization Energy. Besides the electrostatic energy, we have calculated the polarization energy of the solute. It is not necessary to calculate explicitly the polarization energy of the solvent by the solute since, in a continuum model, this is already included in the electrostatic energy, due to the fact that the dielectric constant (by its temperature-independent part) takes into account the polarizability of the solvent. In the same way as for the previous term, we use the "charging-parameter" method!' The variation of the free enthalpy of polarization of the solute AG@ is obtained as the reversible work of charging the polarizability of the solute. In the elementary step where the charging parameter [ vary from f to f df, the atomic polarizabilities a,vary from fa, toJf + df)a,; hence da, = a,df. We assume that the reaction field 6' does not vary during the charging process of a,;hence

+

and

(45)

but

Note that, here, the charging process has introduced no additional factor, due to the assumption above (namely, neglecting is the potential created by the charge distribution p 0 ( 3 at the point hence, finally

(61) Yosom, S. J.; Owens, B. B. J. Chem. Phys. 1963, 39, 221. (62) Constanciel, R. Theor. Chim. Acta 1986, 69, 505.


Langlet et al.

the variation of the reaction field during the charging process for the polarizability). Within our methodology, the reaction field &(Fi) is simply calculated as the electric field created by the charge surface density

where by definition

a:

is the atom-atom distribution function for a given value of r,,”. As in ref 42 the treatment used in the present model for the calculation of dispersion-repulsion energy corresponds to eq 52-58 with p 2 ( i j ) replaced by an approximate distribution function g,,, corresponding approximately to a solute hard sphere with radius R: embedded in a solvent of hard spheres with radius Rr (where radii R: and RY denote the respective van der Waals radii). This was achieved4*by evaluating calibration constants KFv (for dis(for repulsion) from a packing of hard spheres persion) and , with radii RY (simulating the solvent) around a hard sphere RW (simulating the solute) for each ~ , pair. u Furthermore, as in ref 42, we neglect (i) the reorganization of the solvent associated with the switching on of the dispersionrepulsion energy, so that

A@“ defined by eq 45 is nothing but the mean polarizatio? energy of the solute by the solvent. The neglect of variation of &l under the charging process of the atomic polarizabilities a, leads to neglect of the reorganization energy of the solvent due to the polarization of the solute, but we think that this effect is very small. Let us stress that the polarization energy calculated above is an approximation, for two reasons: (i) the reaction field we calculate is not really the exact one but an approximate one; (iil at any time t , the electric reaction field at any point F!, namely, &’(F,,t)is the field created at this point by the solvent molecules occupying their position at time t. The exact polarization energy at time t is given by (v2v)l

= -(1/2)a,((B’(F1,t))’)

(47)

and what we call reaction field is &(F,)

= ($’(?J))

(48)

In eq 45 appears the square of a mean value of the electrostatic reaction field, whereas in eq 47 appears the mean value of the square of the electrostatic field so that eq 45 should really be written as AGP’ = ( Q!‘iv)l = -( 1 / 2 ) ( ~ ,(Z’(7,J)))’ (

+ A&”

(49)

with A&‘’ = ( (Z’(F,,t))’)- ( ( Z’(?,,t)))’

(50)

A&’, the standard deviation of &(i,,t) with regard to its mean value has been neglected in eq 45. Since the polarization energy is weak with regard to the other components of energy, we do not feel, however, that the approximations made are important. (c) Expression of the Dispersion-Repulsion Part of SoluteSolvent Energy. The exact average dispersion-repulsion energy is

( U ) = J . . . J U ( i l ) f(il) dil

(51)

where U(Q)is the total interaction energy, Q stands for the set of all coordinates of the molecules, andflil) is the weight function. U(il) is a sum of pair interactions

gpV(r@U) l:2(qi,q,) dqJ

(58)

e,

A

G

N

~(vf-;,”)l ~

(59)

and (ii) the entropy part of the dispersion-repulsion energy, so that

m-RA N

G

N

(e:)l ~

~

(60)

Approximation i may be valid for polar molecules in which the essential part of the reorganization energy of the solvent proceeds from the electrostatic interaction rather than from the dispersion repulsion energy, but it must be less justified in the case of nonpolar solute and further studies of this point could be necessary. The evaluation of the reorganization energy of the solvent, due to the dispersion-repulsion interaction, appears less advanced than in the case of electrostatic, since there is presently nothing similar to the dielectric model for this interaction. An important point is that, by contrast with the case of the does not vanish (even a electrostatic interaction, imately) for [ = 0. We shall tentatively assume that ( s-sv linear function of (, namely

(e:)(

=

(eiF)o+ t[(@4:)1

J-!?;

-

(uP-;yR)ol

(61)

-

(ei:)oI

(62)

and by inserting this into eq 56, we get

AGDR =

( e F ) o + X[(@iF)i

AGbR = Yz[(%F)o

+ (ei:)11

(63)

(52)

From a practical point of view, presently, the dispersion-repulsion energy is calculated as the sum of roducts of energy integrals by calibration constants Kf: and

where q, (and q,) represents the set of the three positions and three angular variables necessary for specifying the position of i (and j ) . Hence, eq 51 can be written as

(64)

U(Q) = CVD-R(q,,ql) 121

( u, =

( N / v)

1VD-R(q19qJ)P2(q~3qI)dqJ

(53)

where N is the total number of molecules, Vis the volume, and p2(q,,ql)is the distribution function. The integral (53) only depends on the relative position of j with respect to i. When p R ( q , , q Jis) expressed as a sum of atom-atom ( W J ) terms UD-R(qid7,) = C C q ; R ( r , , ” ) &il Y t J

we get

with

(54)

6

where i is a summation on each atom of the solute molecule, j is a summation on each distinct species of atoms of solvent, nl is the number of atoms of type j in the solvent molecule, and V$‘ is the effective volume of the solvent, i.e., the molar volume divided by Avogadro’s number. e t and ef:are respectively the dispersion and repulsion potentials of pairs of atoms i and j ; we have used potential functions of Kitaygorodski.66 Then, with Ostrogradski’s formula, the volume integrals of eq 64 are transformed into surface integrals which are calculated numerically by using Korobov’s grid. For more details of this part of energy see ref 42. In the framework of the present procedure, it must be acknowledged that the obtention of accurate values of the ther(63) Langlet, J.; Claverie, P.; Carbon, F.; Boeuve, J. C. I n r . J . Quantum Chem. 1981, 20, 299. (64) Ben Naim, A. J. Phys. Chem. 1978, 82, 792. Ben Naim, A. Hydrophobic Interactions; Phenum: New York, 1980. (65) Ben Naim, A.; Marcus, Y. J . Chem. Phys. 1984, 81, 2016. (66) Kitaykorodski, A. I . Tetrahedron 1961, 14, 230.



modynamic quantities related with the dispersion-repulsion part of the potential energy represents a nontrivial problem: (1) For obtaining the calibration constants Kf: and appearing in eq 64 above, the procedure used a distribution function pertaining to a mixture of spheres (corresponding to the various atoms) instead of real solvent molecules; thus we cannot pretend or ( pR) However, since ( pR) to reach exactly either ( PR), is obtained as an integral involving the distribution function, we may expect that its value is not very sensitive to the details of this could be reasonably function, so that our evaluation of ( pR) accurate in spite of the approximation. A test of the validity of this expectation consists in evaluating the vaporization energy as half the average interaction energy of one molecule (the solute) with the surrounding ones (the solvent), as is done in section V.A. Indeed, for nonpolar molecules, this evaluation reduces to ( pR) and we therefore have here a convenient criterion for the validity of our approximation. (2) Equation 59 is approximate in so far as it assumes that the distribution function is independent of the charging parameter [. This point was theoretically discussed in ref 42. In fact, the comparison of calculated and experimental values of AG in the present work (section VI(B)) will enable us to check the practical validity of this approximation. 111. Attempt To Take into Account Solvation Sites The macroscopic dielectric constant does not give a perfect representation of the statistical behavior of solvent molecules in the immediate vicinity of the solute molecule. This model should not give the "solvation sites", Le., the sites of the solute molecule that interact very strongly (by H bond for instance) with some solvent molecules. The introduction of a factor f (different for polar and nonpolar solutes) that multiplies the electrostatic free energy has given satisfactory results as concerns the enthalpy of vaporization of several types of molecules,42 but this procedure does not represent solvation sites of the molecule. Therefore, we have searched for a systematic way to define a modified van der Waals radius for each atom of the solute molecule. A possibility for reducing the van der Waals radii is to take into account the fact that the cavity containing the solute is compressed due to the electrostatic interaction between the solute molecule and the solvent. We have tentatively adapted formulas given by Whalley6' for a sphere containing charges, dipoles, and quadrupoles embedded in a dielectric medium. This author used the electrostatic free energy of solvation AG,, given by Kirkwood52 and obtained the electrostatic volume change AVa of the system by differentiating AGsolvwith respect to pressure, P. AV,, thus obtained can be subdivided into two parts: the first one due to terms in dc,/aP, which is a change of dielectric volume (thus its electrostriction), is negligible according to Whalley and the second one is a change of volume of the cavity caused by the force exerted on the cavity by the dielectric and is expressed as

where a In a/aP is -5 X lod bar-' as in ref 64, N is Avogadro's number, and q, p, and 0 represent respectively the charge, dipole, and quadrupole of the sphere of radius a, embedded in a medium with dielectric constant cr. Equation 65 can be written as 1 a h a AV,,(cav) = -N-[AG(O) 3AG(1) 5AG(2)] (66) 2 ap

+

+

where AG(O), AG(l), and AG(2) are the electrostatic free energy due to the charge q ( I = 0), the dipole p ( I = l ) , and the quadrupole 0 ( I = 2) carried by a sphere of radius a,. (67) Whalley, E. J . Chem. Phys: 1963, 38, 1400.

Once AVa has been obtained, it is easy to calculate the modified radius a of the cavity by AV = (4/3)?r(a3 - ao3)

(67)

Thus a = (ao3

+ (3/4~)Av)'/~

(68)

Equation 66 obtained for a single sphere can be easily transferable to our cavity built from the union of van der Waals spheres (see Figure 1). For each atom i, we indeed define contributions AG(i,l) as follows: we start from eq 13 and 39, which give Ace' =

~2~0(s3P,0,(s3 S dS

(69)

Inserting the decomposition of P,,(eq , 40) we obtain the corresponding decomposition AGC1(i,O)+ AGel(i,l)

AGel =

+ AGe'(i,2)

(70)

i

where AGel(i,l) = y210(s3P(i,I,s) ds'

(71)

Now for each atom i, we can use eq 66, with AG( 1) given by AGel(i,l) defined in eq 70. Then the modified van der Waals radius of atom i, namely, R'I is calculated according to eq 67, where a and a. are replaced by R'Y and RY, respectively. We define a modification factor of each van der Waals radius as

ART) = R'T/RT

(72)

For reason of simplicity, in this part of the calculation, we have used only atomic multipoles 1 (1 = 0, 1, 2) obtained from our multipolar d i s t r i b ~ t i o n . ~ ~ Furthermore, as an attempt to represent the hydrogen bonds when present, a multiplicative factor Cihas been introduced in the calculation of R'I (by eq 68):

where q, is the net charge on atom i, N,(i) the number of electrons on the valence shell of the atom i, and P, the electronic population of the valence shell of atom i. The factor P,/Nv(i) has been used in the calculation of the repulsion part of interaction energy between discrete molecules as a way of taking into account the reduction in repulsion for short distances accompanying the decrease in the electronic density. We expect that the factor C, could represent in an analogous way the higher compressibility of a hydrogen bearing a high positive charge. In Table I we give such modified van der Waals radii calculated for some nonpolar molecules (cyclohexane, CC14) and polar molecules (acetone, dimethyl sulfoxide and dimethylformamide), with and without the factor C, of eq 73 (column 1 or 2, respectively). The two results are nearly equivalent. On the contrary, in water or methanol, the modified van der Waals radius of the hydrogen entering the H bond is very different according to whether the factor C, is introduced or not. Hence we recommend introducing this factor in all cases. IV. Test of the Validity of Our Approximate Fictive Charge Density Q for Calculation of Electrostatic Energy It may be suspected that three major elements will influence the validity of our approximations: (i) the presence of some strongly charged atoms in the central molecule, (ii) the size of this central molecule, and (iii) the dielectric constant of the medium outside the cavity containing the solute. To investigate these effects we have chosen five molecules: acetone and dimethyl sulfoxide (same size, nearly same geometry, but different dipole moment; p = 2.88 and 3.84 D, respectively), water, methanol, and 1-propanol (same charged group OH but different size)


of acetone, AG; is slightly overestimated with regard to AC;& for low values of the dielectric constant e, < 7.4, while the two results coincide exactly for e, > 7.4. In the case of DMSO, the difference between the two calculated values remains constant in the whole range of the dielectric constant. In the case of water and methanol, the two calculated values of AGel are very close in medium of very low dielectric constant e, = 2.28 but AG: is slightly overestimated with respect to AQ&, and the difference between the two terms slightly increases with the dielectric constant. In the case of 1-propanol,the behavior of AG; and AGgc is similar to that encountered with acetone. On the whole, the difference between AG$ and AGgc is never larger than 15% (water in water), so we think that the simplified methodology may be used for the study of solvent effects at least in a first approximation, or for very large solutes. Note that when very accurate results are required in the calculation of the electrostatic energy, the approximate evaluation of u might be used as a starting point of the iterative calculation of eq 28. In the next section, we shall see that, in most cases, our calculated vaporization energies are very close to experimental ones. In a later paper we will show that good agreement with experiment can be obtained in the calculation of the effect of nonpolar solvents upon the H-bond formation of the nucleic acid bases. At this point, let us add one more remark concerning the usefulness of going through the relatively elaborate steps of our procedure, namely, the use of a "multisphere" cavity to represent the molecular shape and of a fictive charge density on its surface, rather than a simple spherical cavity of appropriate radius enclosing a simplified charge distribution (due to the fact that such a sphere cannot match the molecular shape, it would then generally not be possible to use the set of multipoles on atoms and middles of bonds used here). The fact that the use of a spherical cavity is not very realistic for complex nonsymmetrical molecules has been discussed many times (see, for example, ref 59 and references cited therein and more recently ref 60 and 62). Indeed, it was explicitly shown60 that the numerical values of the electrostatic energy strongly depend on the radius adopted for the sphere and on the choice of its center. On the other hand, it has been shown63 that, while the molecular charge distribution can be satisfactorily represented by sets of multipoles (charge, dipole, quadrupole, for example) located at seueral centers (atoms and bond centroids for example), it cannot be so, in general, by using a set of mul-

TABLE I: Modification Factor of the van der Waals Radii in Various Molecules' molecule cyclohexane carbon tetrachloride

atoms

C H C CI

c (C=O) c (CH,)

acetone

0 dimethyl sulfoxide

H (CHd S (S=O) C (CHJ 0

dimethylformamide

H ((343) C (C=O) c (CH,)

0 N

€3 (CH,) H (H-C=O) 0

water

H

c (CH,)

methanol

0

H (CHJ H (0-H)

f(RX1) 0.998 0.998 0.994 0.996 0.999 0.990 0.989 0.990 0.992 0.998 0.937 0.979 0.987 0.995 0.995 0.999 0.995 0.999 0.978 0.964 0.994 0.976 0.990 0.976

Langlet et al.

ARX2) 0.998 0.998 0.979 0.993 0.942 0.972 0.965 0.964 0.957 0.976 0.909 0.955 0.943 0.942 0.961 0.969 0.903 0.928 0.923 0.837 0.980 0.941 0.890 0.810

'f(R;U)(l) is calculated according to eq 72 andf(R;")(2) is calculated as C,flR;U)(l) with Ci given by eq 73.

plunged in various dielectric media, with the dielectric constant ranging from 2.28 to 78.0. Then, for each molecule, we have calculated Ace' as a function of the dielectric constant of the outside medium within the present methodology, A@, and within that proposed by Huron et AG&. All calculations have been performed with nonmcdified van der Waals radii and with our atomic monopoles and dipoles only since the method of ref 59 has been developed for such a charge distribution only. The electrostatic free enthalpies A G ' and AGAc are collected in Table IIA. The data show monotonous variation of A E ' and Aff& with the dielectric constant, but the two series of values do not coincide exactly; it is observed in particular that in the case

TABLE 11: Values of Electrostatic Free Enthalpy AG' 'r

molecule acetone

AGHc AG, A AGHc AG, A AGHc AG, A AGHC AG,

2.28"

4.80b -1.56 -1.71 0.15 -3.43 -3.71 0.28 -3.38 -3.06 -0.32 -2.27 -2.08

propanol

AGHc AGc A

-1.04 -1.24 0.20 -2.28 -2.56 0.28 -2.20 -2.15 -0.05 -1.47 -1.45 -0.02 -1.24 -1.04 0.20

water

AGs A AGS A AGs A

-0.64 1.56 -0.22 1.25 -0.07 1.17

dimethyl sulfoxide water methanol

A

methanol propanol; R, = 3.35

A

7.4OC

13.10"

20.0oC

36.73f

46.489

78.00''

-1.64 -1.80 0.16

(A)' -1.74 -1.86 0.12 -3.84 -4.10 0.26 -3.79 -3.36 -0.43 -2.56 -2.29 -0.27 -1.86 -1.99 0.13

-1.90 -1.97 0.07 -4.18 -4.42 0.24 -4.15 -3.60 -0.55 -2.80 -2.46 -0.34 -2.05 -2.14 0.09

-1.97 -2.02 0.04 -4.35 -4.57 0.22 -4.32 -3.72 -0.60 -2.92 -2.54 -0.38 -2.15 -2.21 0.06

-2.03 -2.06 0.03 -4.48 -4.69 0.21 -4.45 -3.81 -0.64 -3.02 -2.61 -0.41 -2.22 -2.27 0.05

-2.04 -2.08 0.04 -4.52 -4.72 0.20 -4.49 -3.83 -0.66 -3.04 -2.62 -0.42 -2.24 -2.28 0.04

-2.07 -2.09 0.02 -4.57 -4.71 0.20 -4.55 -3.87 -0.68 -3.08 -2.65 -0.43 -2.27 -2.3 1 0.04

-1.02 2.36 -0.35 1.92 -0.10 1.54

(BY -1.15 2.64 -0.39 2.17 -0.12 1.74

-1.27 2.88 -0.43 2.31 -0.13 1.92

-1.32 3.00 -0.45 2.47 -0.13 2.02

-1.37 3.08 -0.46 2.56 -0.14 2.08

-0.47 2.57 -0.14 2.14

-1.40 3.15 -0.47 2.61 -0.14 2.17

-0.19

"CCl4. bCHClp. 'THF. "Pyridine. CAcetone. fDMF. gDMSO. *H20. 'Ace' calculated within the methodology of ref 59 (AGeAc) and with the method presented here (AE'). A represents the difference APAc - AG?. 'AG$ calculated with the exact formulas for the surface charge u but assuming the cavity shape as a single sphere corresponding to the experimental molecular volume. A represents the difference ACic - AG:'. All values in kcal/mol.


Improvements of the Continuum Model tipoles (even up to higher orders) located at a single center. As an illustration of the situation, we have calculated the electrostatic free enthalpy of water, methanol, and propanol represented by spheres with radii corresponding to the respective molecular volumes used in the calculations through our own method; at the center of the sphere we placed a dipole obtained from the charge distribution used in these previous calculations. The results from the %ngle-sphere” calculations are given in Table IIb. They clearly appear as severely underestimated (by 66% for water in water, for example) when the sphere corresponds to the experimental molecular volume. Unrealistic reduced values of the cavity radius would have to be adopted for getting reasonable energy values (and these values should probably vary, for a given solute, according to the solvent). It therefore appears more realistic to use the “multisphere” van der Waals cavity. This procedure seems accurate enough provided that the van der Waals radii are suitably adjusted. In fact, this could compensate for the part that the solvent is approximated as a continuous dielectric. One might of course conceive of further improvements (e.g., introducing a variable dielectric constant in the vicinity of the solute). V. Expression of Thermodynamic Quantities of Solvation In the calculation of the thermodynamic quantities of solvation AGmI,and AHmlvr we have taken into account the variation of the solvent volume A P l on the cavitation step (see eq 3) but not the negligible variation of solvent volume A P 2 due to the electrostatic pressure, so these thermodynamic quantities can be directly compared to the ones given by Ben Naim64 Le., AG*so~vand AP&. This author has defined the process of solvation of a solute s in a solvent 1 as the process of transferring the solute s from a fixed position in an ideal gaseous phase (8) into afixed position in the liquid (l), with the process being carried out at constant temperature and pressure, all other components of the system remaining also unchanged. Though the procedure is not an experimentally feasible process, Ben Naim has developed a simple statistical mechanical argument which provides a connection between the coupling work of the solute s to the liquid 1 w(s/l) (which is nothing but the insertion free energy AG*mlv)and experimental quantities. Then the other thermodynamic quantities of solvation have been obtained by standard procedures. It has to be noted that the starred thermodynamic quantities introduced by Ben Naim do not contain terms due to the compression of the gas from the gas phase to the liquid phase. In particular, AH*mlvdoes not contain the classical term P( V , - V,) N RT; AH*,], only contains the term P ( A P , + A P , ) defined above. Therefore, AH*solvis very close to the solvation internal energy AEsolv,which in pure liquids is obtained from the heat of vaporization obtained at constant V. In this section we briefly recall and discuss the formulas giving within our methodology the various thermodynamic quantities of vaporization in the case of a pure liquid: (A) enthalpy of (B) free enthalpy of vaporization (AG*,,,), vaporization (APvap), and (C) entropy of vaporization (AS*vap). ( A ) Enthalpy of Vaporization. The solvation process of vapors condensed into their liquids at a given temperature T is just the opposite of the vaporization process. In the special case of pure liquid, the solvation energy may be interpreted in terms of the average binding energy of the solute s to its surrounding plus the polarization energy of the solute by its surrounding: -AEvap = AEsolv =

1/2(

us,,) 1 + ( W 1)

(74)

where (U,,,) I is the ensemble averaged interaction energy of a single molecule s with all the other molecules of the solvent, namely We can also calculate APEmI,(or AHsolv),in the framework of the solvation process: creation of the cavity and switching on of the different parts of energy (dispersion, repulsion, electrostatic, and polarization). We have shown in section I1 that our calculated values represent solvation enthalpy rather than solvation internal

A&”

creation of cavity

II

Introduction of the solute into the cavity

IV

F,=1

Figure 2. Scheme of the solvation process in nonpolar liquids: AH*solv = AH*,, + Hi, + H11.

energy. Thus using Ben Naim’s insertion quantities62AH*, AS*, and AG* which leave out the term PAY = R T , we obtain

+A

-AH*,,, = AH,,,

P R + AHe‘

+ AH@

(76)

Then, within our approximation, by use of 59, 60, and 15, eq 76 may be written as -AH*”,, = AH,,,

+)!A!@(

1

+ (&,)

1

+ AU&, + AIP” (77)

Note that this equation involves the same approximations as eq 59 and 60, namely, the neglect of the reorganization of the solvent associated with switching on of the dispersion-repulsion potential. Combining eq 75 and 77, we obtain (since AH*,,, N AEvap) AH*,,, = AHav AG:,, AHPOI - 2AGP”I (78)

+

+

In the case of nonpolar (or weakly polar) liquids, the electrostatic part of the solvation energy is negligible, and the same is true of the reorganization of the solvent due to the switching on of the electrostatic part of energy. Thus eq 75,76, and 78 reduce respectively to -AH*,,,

N

-AEvap =

+A

-AH*,,, = AH,,

Yz(

+ (@!?$),

P R = AH,,

AH*vap =

(79)

1

(80) (81)

m c a v

In fact, when we consider the correct solvation process as described in Figure 2, we see that eq 80 should actually be written as -AH*,,, = AHcav +

[(%%)I

- (%Z)OI

+

(82) (e!)]

or -AH*vap =

w c a v

+ A%&

+

(@!
I

(83)

Equation 83 differs from 80 by the introduction of the solvent reorganization energy A&.. during the switching on of the dispersion-repulsion energy. With this correction, eq 8 1 becomes AH*vap= AH,,,

+ AV;:

(84)

The study of nonpolar (or weakly polar) liquids in comparison to experimental values can thus provide a test on the accuracy

Langlet et al.


of our calculated average dispersion-repulsion energy through eq 79. On the other hand, two alternative tests are offered by using eq 81 or its nonapproximate version (84): should the dispersion-repulsion reorganization term be really negligible, eq 81 could serve to test the calculated cavitation enthalpies. Conversely, and more likely, if this term is not negligible, a reasonably correct evaluation of the cavitation enthalpies and use of experimental vaporization energies in (84) could provide an evaluation of Ai%&. Concerning the cavitation enthalpies, we may expect a reasonable accuracy from the formulas given by Sinanoglu (eq 2a,b) since the parameter Kf(b/a) which fits the planar surface energy to highly curved microscopic dimensions has been obtained from experimental solubility and the mixing data of asymmetric syst e m ~ As . ~ concerns ~ Pierotti’s formulas (eq 3-6), the ability to determine from experimental data the solubility of a hard sphere in a real solvent has allowed checking of the adequacy of the scaled particle theory44to calculate AG-,, and a strong argumentation ~~ on the validity of eq 3 and 4 has been given by P i e r ~ t t i .Furthermore, it has been well established that the solubility of gases in organic unassociated solvents as well as in water was adequately treated by his formulas. At this stage, we therefore prefer the second alternative, namely, considering as reasonable our theoretical evaluation of AH-,and extracting from eq 84 an estimation of A=& with the purpose of testing the theoretical evaluation of this term in the future. ( E ) Free Enthalpy of Vaporization. Within our model AC*vap is calculated according to -AG*vap = AG*solv = AGcav

+ AGC1+ AGP”’ + AGD-R

(85)

As was done for AH*,62 we have used an asterisk (AG*) to denote that the calculated free enthalpy does not contain the term PAVdue to the compression from the gas phase to the liquid one.62 Let us reiterate that AGDR has been calculated according to eq 59 and not according to the more exact eq 63. Since this approximation consists in writing

a comparison between calculated and experimental values of AG*,, will allow us to give an estimation of (et),, (see the next section). (C) Entropy of Vaporization. Within our methodology AS*,,, is calculated by

-AS*vap= AS*,Iv = AS,,,

+ ASe’ + ASP”’

(87)

As in ref 42 and 46, we have neglected the dispersion-repulsion component of the entropy. Comparison between the calculated values and experimental ones6’ has allowed us to check this approximation. In the same way as Ben Naim, the ”insertion” entropy AS*vap is given by -AS*vap = AS*,,,, =

a(AG*,lv) aT

Thus it contains neither the change of entropy due to the compression of the gas from p5 to p\ nor the terms including isobaric thermal expansion of the gas and liquid phases, namely, a% and ab. The complete vaporization entropy would be given by -ASvap

=

ASsolv

=

AS*solv

+ R In ( P % / P L )+~ R(ah - ak) (89)

VI. Results and Discussion for Nonpolar and Nonassociated Polar Liquids The method described above has been applied to the study of vaporization energy for a few nonpolar and nonassociated polar liquids at T = 288 K. We shall discuss first the cavitation thermodynamic quantities calculated by the three different methods of section 11, and then we shall consider the calculated

TABLE III: Thermodynamic Quantities of Cavitation Calculated in Terms of Macroscopic Surface Tension 7 (Lines 1)’ Following Sinanoglu’s Theory (Lines 2)’ and Using Pierotti’s Formulas (Lines 3)”

cyclohexane benzene toluene carbon tetrachloride chloroform aniline nitromethane nitrobenzene acetone dimethylformamide dimethyl sulfoxide

5.89 4.78 7.61 5.94 4.9 1 7.57 6.82 5.65 8.34 6.05 4.90 7.84 5.43 4.20 7.47 9.54 8.45 11.52 5.92 5.60 7.00 10.61 7.64 10.20 4.27 3.80 6.43 7.01 9.20 8.05 8.59

24.29 11.30 -1.27 24.42 11.40 -0.33 22.32 11.74 -2.14 24.26 11.57 -1.07 17.92 8.20 -1.14 18.32 20.56 -5.39 21.61 7.00 -2.80 25.87 11.71 -4.37 20.20 10.73 0.27 11.47 -4.46 20.83 -6.90

13.13 8.15 7.23 13.22 8.31 7.46 13.47 9.15 7.70 13.28 8.35 7.52 10.77 6.67 7.13 15.00 14.58 9.9 1 12.66 7.69 6.16 18.31 11.13 8.90 10.29 7.00 6.51 10.43 (7.87) 14.26 6.53

AGav and AHavare given in kcal/mol. All values are calculated at T = 298.15 K. ASs,,, is given in eu (cal/(mol.deg)).

values of the quantities AH*vap, AG*vap,and AS*vap. ( A ) Cauitation Thermodynamic Quantities. These quantities are reported in Table 111. For each liquid, lines 1, 2, and 3 correspond to values calculated according to eq la,b, 2a,d, and 3-6, respectively. We have limited ourselves to liquids for which the necessary parameters are available in the literature (K:( 1) and K; in eq 2a-d and a in eq 3 and 4. Table I11 shows that (a) The values of AH,,, calculated by following Sinanoglu’s formulas (line 2) and Pierotti’s formulas (line 3) are relatively close, while the values calculated in terms of the macroscopic surface tension (line 1) is appreciably higher. Clearly, the transformation of the macroscopic surface tension y to microscopic dimensions is of great importance. (b) AGav calculated in terms of the surface tension (with or without the corrective factor) is everywhere smaller by up to 20% than that calculated by using Pierotti’s formulas. We have emphasized, in section II.A, that the essential difference between Sinanoglu’s and Pierotti’s formulas is the presence in the latter of the constant term KO.Indeed we observe that the differences obtained are approximately constant positive and rather important so that the differences obtained between AGav calculated according to eq 2 and 3 should essentially disappear if the term KO is introduced in eq 2c. (c) AS,,, is large and positive when calculated in terms of surface tension (with and without the correction factors) while it is small and negative when Pierotti’s formulas are used. We think that the values of AS,,, obtained in term of the surface tensions are not compatible with the negative values of AS*vap obtained by Ben Naim6* (vide infra). Thus in the next sections we shall only use the thermodynamic values of cavitation calculated with Pierotti’s formulas. ( E ) Vaporization Enthalpies, AH*,,p. (a) Table IV gives the vaporization energy calculated according to eq 75. For each liquid, we have given two results: the first obtained with nonmodified van der Waals radii and the second one with van der Waals radii

The Journal of Physical Chemistry, Vol. 92, No. 6,1988


cyclohexane (2.0)b benzene (2.28)b toluene (2.37)b carbon tetrachloride (2.23)b chloroform (4.80)b dimethyl ether (5.02)b diethyl ether (4.23)b tetrahydrofuran (7.38)b aniline (6.77)b pyridine (13.10)b acetone (20.70)b nitromethane (36.56)b nitrobenzene (34.82)b dimethylformamide (36.73)b dimethyl sulfoxide (46.48)b

-7.94 -7.94 -7.90 -8.05 -8.78 -9.60 -1 1.40 -11.40 -10.05 -10.13 -3.63 -3.76 -5.94 -6.17 -7.50 -8.38 -1 1.93 -12.99 -8.39 -9.16 -5.33 -5.78 -5.73 -6.18 -12.22 -12.69 -7.90 -9.10 -8.60 -9.43

0.46 0.46 1.07 1.10 0.98 1.38 4.14 4.14 3.79 3.82 0.21 0.25 0.29 0.36 0.91 1.47 2.57 3.52 1.48 2.03 0.54 0.77 1.34 1.80 2.99 3.43 1.23 2.09 1.86 2.74

-0.08 -0.08 -0.49 -0.51 -0.35 -0.43 -0.14 -0.15 -0.88 -0.93 -1 .oo -1.10 -0.60 -0.67 -1.56 -1.67 -1.49 -2.04 -1.65 -1.96 -2.21 -2.46 -3.85 -4.21 -3.16 -3.45 -2.38 -2.64 -3.96 -4.82

0.00 0.00 0.00 -0.01

0.00 -0.01 0.00 0.00 -0.05 -0.09 -0.07 -0.08 -0.04 -0.05 -0.14 -0.15 -0.07 -0.12 -0.19 -0.24 -0.07 -0.08 -0.19 -0.22 -0.21 -0.29 -0.10 -0.10 -0.54 -0.76

-7.56 -7.32 -7.47 -8.15 -8.66 -7.40 -7.41 -7.12 -7.33 -4.49 -4.69 -6.29 -6.53 -8.29 -8.73 -10.92 -1 1.63 -8.75 -9.33 -7.07 -7.55 -8.43 -8.81 -12.60 -13.00 -9.15 -9.75 -11.24 -12.27

-7.50

1627

-7.78 -8.75

-7.23

-7.43 -7.06

-3.79 -5.92

-6.45

-6.70

-6.79

-12.80 -10.04 -7.03

-7.30

-8.85 -12.14 -10.62 -1 2.05

and AH*,,, are experimental values obtained respectively from the heat of vaporization and from Ben Naim's process. The table also gives the different parts of the solute-solventenergy ( with X = D (dispersion),R (repulsion), and el (electrostatic), and (TIrepresents ) the solute polarization energy. For each molecular species, this table gives two results: the first one (line 1) has been obtained with nonmodified van der Waals radii while the second one is obtained with van der Waals radii modified according to eq 68. All values are in kcal/mol. bValuesin parentheses are dielectric constants (e,). TABLE V Vaporization Enthalpies A*.,(l)

(A) and AH*,w(2) (B) Calculated According to (A) Eq 76 and (B) Eq 78, Respectively'

AH&( 1) calcd by eq 76 cyclohexane benzene toluene carbon tetrachloride chloroform aniline acetone nitromethane nitrobenzene dimethylformamide dimethyl sulfoxide

eR m1

AHPOI

p'

-M*vap(l)

-14.95 -13.88 -16.44 -14.53 -12.60 -18.94 -10.04 -8.77 -18.52 -14.00 -13.39

0.00 -0.03 -0.02 0.00 -0.20 -0.21 -0.13 -0.32 -0.32 -0.15 -1.13

7.23 7.46 7.70 7.52 7.13 9.91 6.51 6.16 8.90 7.87 6.53

-7.87 -7.29 -9.48 -7.35 -7.38 -12.28 -7.32 -8.65 -14.52 -9.79 -14.27

-0.15 -0.84 0.72 -0.34 -1.7 1 -3.04 -3.66 -5.72 -4.58 -3.51 -6.28

ACb, 0.00 0.18 0.14 0.00 0.15 1.04 1.26 2.70 2.32 1.77 3.36

AH*,,,(2) calcd by eq 78 A(po1) AH*v0p(2)

0.00 0.02 0.00 -0.05 -0.02 0.03 0.03 0.12 0.26 0.05 0.39

7.23 7.66 7.84 7.47 7.26 10.98 7.80 8.98 .11.48 9.69 10.28

AHSBN 7.59 7.78 8.75 7.43 7.06 (12.80) 7.30 (8.85) (12.14) (10.62) (1 2.05)

"AHx (with X = D-R, el, pol, cav) represents the different parts of the total vaporization enthalpy: dispersion-repulsion, electrostatic, polarization, and cavitation. AV$,v represents the reorganization energy of the solvent due to the switching on of the electrostatic solute-solvent energy. A(po1) represents the term AHPOI - ZAG@ in eq 78. All values in kcal/mol. In the last column, the values in parentheses are heats of vaporization.

modified according to eq 68. We observe that for liquids with very low dielectric constant er the results are very close, but when e, increases, the modification of van der Waals radii leads to an increase (in absolute value) of the solutesolvent mean electrostatic energy ( e s v , 1. On the whole (except for ether molecules) our results agree very well with experiment. The slight overestimation obtained for ether molecules can proceed either from an overestimation of the solutesolvent electrostatic interaction energy via the approximation made upon the surface density a or from an overestimation of the solutesolvent dispersion interaction energy: it is possible that the parameters used for an ether oxygen are not the best ones (the same parameters have been used for ether, hydroxylic, and ketonic oxygens). The good agreement obtained between calculated and experimental values of AHvapfor liquids with low dielectric constant shows the accuracy of the mean solutesolvent dispersionrepulsion energy obtained within our methodology. This accuracy certainly proceeds from the introduction of appropriate calibration

Kf:

constants $ and (see ref 42 for more details). (b) Parts A and B of Table V give vaporization enthalpies calculated according to eq 76 and 78 respectively (which neglect the reorganization of the solvent due to the switching on of the dispersion-repulsion potential). They have been obtained with van der Waals radii modified according to eq 68, and AH,,, is calculated according to Pierrotti's formulas. Our results indicate the following features: (i) A P ' increases when e, increases. This term, very small in cyclohexane (-0.15 kcal/mol) reaches -6.28 kcal/mol in DMSO. (ii) AH@ is always small (-1.13 kcal/mol in DMSO). (iii) AC$-,, (reorganization energy of the solvent due to the switching on of the electrostatic energy) is negligible for 2 S er d 4.80; then it increases with tr, reaching 3.36 kcal/mol in DMSO. Overall, when M,, is calculated according to eq 78, its major part proceeds from the cavitation enthalpy. As concerns the values themselves of AH*,,, calculated according to eq 76 (Table VA) or eq 78 (Table VB), the comparison


Langlet et al.

TABLE VI: Vnporization Free Enthalpies and Entropy"

-AG*,,, (calcd)

-AG*vap

(exptl)

ASwv

-7.42 -6.83 -8.54 -6.84 -6.15

-4.42 -4.55 -5.17 -4.40 -4.19 -4.20

-1.27 -0.33 -2.14 -1.07 -1.14 0.27 -5.39 -2.80 -4.37 -4.46 -6.90

cyclohexane benzene

toluene carbon tetrachloride

chloroform acetone

-6.15 -9.58 -6.21 -12.06 -7.55 -10.37

aniline

nitromethane nitrobenzene

dimethylformamide dimethyl sulfoxide

AF' -0.23 -1.11 -0.97 -0.64 -2.62 -4.02 -3.35 -5.06 -3.79 -2.91 -4.90

-U*vap

(calcd)

-m*vap

&Pi

(exptl)

ASP"

ASP"

0.00 -0.07 -0.03 0.00 -0.37 -0.17 -0.30 -0.33 -0.10 -0.17 -1.24

-1.50 -1.51 -3.14 1.71 -4.12 -3.91 -9.04 -8.19 -8.25 -7.54 -13.04

-10.61 -10.80 -12.02 -10.16 -9.63 -10.42

-9.11 -9.29 -8.88 -8.45 -5.51 -6.50

-10.02 -3.31 -11.02 -9.74 -8.45 -6.73

" AG*,,,(calcd) and AG*,,&exptl) are respectively the calculated and experimental free enthalpies of vaporization. All values are in kcal/mol. AS-,, AF',and Pi are the cavitation, electrostatic, and polarization parts of the vaporization entropy; AS*,,,(calcd) was calculated according to eq 87. ASP" and ASP" are the dispersion-repulsion entropies estimated according to eq 95 and 100, respectively. All entropy values are in cal/ (mol-deg). with the values obtained by eq 75 (Table VI) indicates that (i) For weakly polar (or nonpolar) liquids (2 < t, d 4.8) except toluene, there is a rather good agreement between these three values. In toluene we notice an overestimation (0.82 kcal/mol) for -AH*,,,( 1) and an underestimation (0.82 kcal/mol) for Ld We . think AH*,,,(2) with respect to the absolute value of G that the overestimation of - A P V a p (1) is due to an overestimation of A P R which, calculated by eq 60, neglects the dispersionrepulsion reorganization energy of the solvent. The same neglect leads to an underestimation of AH*,,,(2). of A=: (ii) When 6 , increases, the differences between AZZT*~,,(~), AH*,,,(2), and AG? : depends on the molecular species: nearly zero for acetone, nitromethane, and dimethylformamide, it reaches 0.65, 1.52, and 2.00 kcal/mol in aniline, nitrobenzene, and dimethyl sulfoxide, respectively. Even in these cases, the neglect of A F R in our calculations seems permissible in an approximate evaluation of AH*,,. (C) Vaporization Free Enthalpy, AG*,,. The quantities reported in Table VI have been obtained under the same conditions (modified van der Waals radii; Pierrotti's formulas for cavitation) as those of Table V. We observe that the calculated values appear overestimated. Let us focus attention on nonpolar liquids with only dispersion-repulsion and cavitation contributions to AG*vap. If we accept that Pierrotti's formulas give reasonably accurate values of AG,,,, we may assert that the overestimation of AG*,,, proceeds from an overestimation of the calculated value of AGDR. should not be Thus the assumption of eq 86 is not valid; ( but rather smaller. An estimate at this stage equal to of ( can be deduced as follows: let us denote AG*(calcd) and AG*(exptl) the calculated and experimental vaporization free enthalpies and and AGZ! the dispersion-repulsion part of the vaporization free enthalpy calculated according to eq 59 and 63, respectively. The vaporization enthalpy has been calculated accordin to eq 85, with the approximate value of AGDR, namely, AGapprox.The corresponding exact value of AG&ft can be written as

(e:)]

e:),

d

AGg$t = AG*,,,(exptl) - AG,,, - Ace' - AGW'

(90)

In the cases where AG*,,,(exptl) is known (cyclohexane, benzene, carbon tetrachloride, chloroform, and acetone), we observe that there is an empirical relation between AGgFo, and AGZ: AG~:,

N

0.8~cg;~~

(91)

From this relation and using eq 63 and 59 we deduce 1/2[(=:)0

+ (=:),I

= O.8(=:),

(92)

leading to

(U?;:)O = 0.6(vfls,R)1

(93)

This evaluation indicates that, within the limits set by the small number of compounds considered, the term ( is appreciably

a:),

e+R)

smaller than ( I and that an adequate theoretical evaluation should be considered in the future. (0) VaporizationEntropy, AS*,p.Table VI gives the different parts of the total entropy calculated according to eq 87. It is observed that (i) In most cases, ASmv represents the major part of the entropy thus calculated. (ii) The electrostatic part of the vaporization entropy is small for tr < 4.8 and increases with t, up to -4.9 cal/(mol.deg) for e, = 46.48. (iii) The polarization part of the entropy is always very small. (iv) The total entropy calculated according to eq 87 is clearly underestimated. This may proceed from two causes: either from an underestimation of the cavitation entropy (which is not very likely since both AG,,, and AHmvcalculated according to Pierrotti's formulas seem to be quite adequate) or from the neglect of dispersion-repulsion entropy. Let us call AS*,,,(calcd) the vaporization entropy calculated according to eq 87 and AS*,(exptl) the experimental value. This last quantity can be written as -AS*,,,(exptl)

= AS*,,,v(exptl) = ASav AS='+ AS@

+

+ ASD-R

(94)

whereas our calculated value is obtained by eq 87. Then an estimation of ASbR can be obtained as = AS*,,,(calcd) - AS*,,(exptl)

(95)

On the other hand, we may try to get a theoretical expression of ASD-R as follows using the relation T ~ D - R= A P - R - A G D R (96) inserting in eq 96, the exact equation for A P R

WR = A@:;

+ (et),1

(97)

and eq 63 giving AGDR, we finally obtain -TASEiR = TU?%: =

+!h[(%?)o

+ (=,")I]

+ A%$.

+ (%:)I

(98)

From this equation, we can obtain a second estimation of ASpR, using the approximate evaluation of (eq 93), and yielding

(e:),

T A S f R=

AV:; + 0.2( vp;,")l

(99)

According to our general standpoint, we can neglect (as in previous sections) A F " and MyR. We have reported in Table VI the values of ASsR calculated by both estimations ASPRand ASPR. Inasmuch as (i) A F R is really negligible and (ii) the approximation (93) is accurate, the values of ASfR and ASpR should be similar. Indeed, the computed values are very close for cyclohexane, benzene, carbon tetrachloride and acetone; the difference that is observed for toluene and chloroform (2.1 and 2.8 cal/(mol-deg), respectively) is, however, relatively small. It appears that the two series of estimated values for ASDR give a good hint for further developments.



VII. Conclusion In the light of the results of the calculations reported above, let us reexamine the various approximations involved in our model and their consequences. (1) As concerns the electrostatic energy, we were faced with three alternatives: (i) to approximate the shape of the cavity through a very simple analytical one (e.g., a sphere) and then solve explicitly the corresponding very simplified electrostatic problem; (ii) to choose an accurate shape of the molecular cavity (obtained as a set of intersecting spheres in our method) and then solve exactly the electrostatic problem by rather heavy computational techniques;59@(iii) to choose the same accurate shape of the cavity while trying to implement a simplified solution of the corresponding electrostatic problem (by using a suitably simplified expression of the fictitious charge density located on the cavity surface which generates the reaction field due to the polarization of the solvent by the solute). Having retained the third alternative, we have shown that, when compared with the exact (numerical) s o l ~ t i o nour , ~ approximate ~~~ value of the solute-solvent electrostatic energy seems accurate enough for practical purposes. We showed furthermore that the details of the shape of the solute cavity are definitely of importance concluding, as Miertus et a1.,60 that calculations concerning large solute molecules embedded in a single sphere lead to rather crude results. (2) As concerns the dispersion-repulsion part of the interaction energy, several questions were raised in section II.c, namely, (i) How accurate is the average dispersion-repulsion energy evaluated (through the fitting of the calibration constants and Kf)by using a solvent distribution function corresponding to a mixture of atomic spheres? (ii) To what extent is the solvent distribution function independent of the dispersion-repulsion cou ling pal ? (iii) rameter; namely, can we write APR = AP-R = ( To what extent can we consider the reorganization energy of the solvent as negligible upon switching on the dispersion-repulsion interaction energy PR? The comparison for several thermodynamic quantities of our calculated values with their experimental counterparts has led to the following respective conclusions: (i) The detailed behavior of the solvent distribution function does not seem essential for the evaluation of the average dispersion-repulsion energy. Actually, in the case of nonpolar liquids (where the dispersion-repulsion part is dominant), we found a satisfactory agreement between the experimental vaporization energies and our values calculated according to eq 79 (which corresponds to describe the solvent through the distribution function corresponding to a mixture of hard spheres surrounding the solute cavity). (ii) By contrast, the overestimated values that we obtained for and are not equal, indicating AG*,,, suggest that that the solvent distribution function should actually exhibit some nonnegligible dependence upon the (dispersion-repulsion) coupling parameter. We therefore proposed a simple rule for evaluating ( from ( based upon the comparison between the calculated and experimental values of AG*vap. As a consequence of this study, we also observed that the entropic contribution associated with the switching on of the dispersion-repulsion interaction is not negligible, and we also proposed a simple recipe for estimating this entropic contribution. (iii) The fact that our calculated vaporization enthalpies are in good agreement with the experimental values (see (i)) suggests that the reorganization energy (associated with the switching on of the dispersion-repulsion energy) should be essentially negligible [despite the fact that the contribution of this reorganization to the entropy is not so (see (ii))]. On the whole, the values calculated for the dispersion-repulsion interaction seem closer to than to ( thus, according to eq 80, we obtain AHER rather than AGDR. A final remark concerns the evaluation of the cavitation contributions to enthalpy and entropy: it must not be forgotten that Pierotti’s formulas (45 and 46) are based upon the scaled-particle theory, which considers hard-sphere interactions only. Thus,

a!)

(e:), (e:)

significant improvements might probably be obtained by taking into account the genuine solventsolvent potential (especially the electrostatic part in the case of polar solvents). This is an important direction now deserving to be thoroughly investigated. Even a t this stage, however, it appears that our procedure may be of practical use in a number of problems where solvent effects must be considered. Appendix I. Calculation of the Reaction Potential V’ from a Fictitious Charge Distribution on the Surface of the Cavity The derivation of the general equations leading (1) to the expression of V’from a surface charge distribution u and (2) to the evaluation of u, is the following: The potential V’ created by a volume dipole 8 is 1 -1 V‘ = --+grad 47rto r dV V’ =

x(diie

1 47r€o

r

- div 8) d V r

(1.2)

where Z is the unit vector normal to the surface. A fictitious surface charge distribution u and a fictitious volume charge distribution p are defined as -Q

= J.8

(1.4)

p = -divp’ (1.5) For a homogeneous dielectric medium (e independent of 7) that does not contain charge distribution, one gets

-

AV= 0 Hence div

8 = div (-grad

(1.6)

V) = -AV = 0

(1.7)

- to) div 2

(1.8)

Thus p

1 = -div p’ = --(e 47r

=0

Hence p = o

Equation 1.3 becomes V’=

ds

(1.10)

The reaction potential V’defined by eq 10 satisfies the condition AV’ = (1 eo)^ = 0 (1.1 1) and the boundary conditions are

e:)1,

(1.12) with

IE

and

Ir

indicating outside and inside the cavity.

In order to calculate u, we write boundary conditions for the total potential V. (1.13) Since the discontinuity of the normal derivative of V only proceeds from the discontinuity of V’ ( Vohas continuous derivatives) eq I. 12 gives

(e,?)

(1.14) Thus combining eq 1.13 and 1.14, one obtains --Q

€0

(1.15)


Figure 3. Definition of the different geometrical elements used for the calculation of the reaction potential.

Hence u is given through the calculation of the normal component of the total electrostatic field inside the cavity. To do that, one has to consider (see Figure 3) (i) a surface element ds'on the surface of the cavity associated with the point M; (ii) a surface element ds and a vector 7 = 3 - 3'as the distance between elements ds and ds'; (iii) a point P defined very close to ds, but inside the cavity as shown in Figure 3. At point P one has to consider (i) the electrostatic field $(a ds') created by all the surtace charge distribution u ds'and (ii) a local electrostatic field &',(a ds) created toward the inside of the cavity by the charge surface distribution of ds:

Langlet et al. €1

+ €2

€1

- €2

cos 8

ds' (1.23)

Appendix II. Relation between the Fictitious Charge Surface Density u and the Normal Component of the Electrostatic Field Created by Multipoles of I inside a Sphere of Radius a We consider a sphere of radius a having a dielectric permittivity to inside which are situated M discrete point charges e, ... em. A polar coordinate system with origin at the center of the sphere may be employed to describe the configuration of the charges e , ...e,. The position of each charge ek is given by three coordinates rk, 8k, &. The electric potential at any point (r, 8, 4) inside the sphere is given by (11.1)

V=vO+V'

with V' as the potential created by the whole charge distribution, which may be developed as a multipolar potential, and V' as the reaction potential given by

V' =

-

+1

~ ~ , r ' ~ y ( c@elm/ os

(11.2)

/=O m--1

p;l(cos %)are the associated Legendre functions and B,, are constants determined by the proper boundary conditions. We write

(1.16)

Y7(8,4) = P;t(cos %)e"'

(11.3)

The normal component at the surface of the sphere of the reaction field is

where

&,, = a v ' / a n

-

+1

b/, = 1=0 c m=-1 c

The reaction electric field on the surface S proper is given by

(11.4)

B/mIa/-1qy%,4)

(11.5)

The nonpal component of the electric field created by the charge distribution inside the sphere is $on

=-av/an

(11.6)

and U

ZrI(uds) = --jj 2CO

(1.18)

The normal component at the surface of the electrostatic field is

denotes the electric field created by the surface element ds with charge density u in its close neighborhood. Using

a 3, XI, avD

=

+

av' Zl,

(1.19)

When the treatment of Kirkwood is restricted to our case, two media, solute neutral (with dielectric permittivity e) dielectric instead of solute neutral dielectric dielectric with ionic strength, the constants Elmand BImare related by

+

and inserting eq 1.16-1.19 into eq 1.15, we get

+

+

Elm

Bh =

(1.20)

(1

+ 1)(€0- e)

toa2" (I + l ) € + l€o

(11.9)

Combining eq 11.8 and 11.10 giving the expression of u

Thus

(11.10) When dealing with

t,

= e/eo, we get

-1_+ % U 1 - e,

€0

one gets cos %

ds' (1.22)

Using the notations of Durand5'

c- c

(11.1 1)

+I

eq 22 becomes

=

/=0 m=-1

u/myiYb$)

(11.12)

J. Phys. Chem. 1988,92, 1631-1639 Using eq 11.9 relying on BI,,,on El,,,, we write

Thus (i) dividing both numerator and denominator by eo and (ii) recalling that the dielectric constant is nothing but t, = €/to, we get (21 l)(€, - 1) = (11.20) (1 l)Er 1

+ +

The expression between brackets in eq 11.1 1 becomes

- €0)

( 1 + 1)E/m - ( 1 +

(I

c0a1+2

+ l)e + le0

1Elm -

eoa/+2

(11.14)

which may be written

I+ 1

&e - €0)

1

According to eq 11.15, the fictive charge distribution 11.12 is

1631

(11.15) ulmof

eq

+

As recalled in section 1I.B the factorfll) calculated by eq 11.19 is different from the corrective factor fMsT of Miertus et Effectively for a dipole 2,for instance, the charge surface density u is exactly given by u

3€, 1 €, - 1 211 cos 8 = 2€, 1 4T a3

7-)

-( +

(11.21)

In the case of dipoles (1 = 1) the corrective factor fMST(1) of Miertus et aL60 is (11.22) Whereas ours is 3(€, - 1) 2€, 1

f(1) = --

+

Hence according to eq 11.8 Ulm

1 (21 + 1)(€ - €0)=Gon(ltm) 47r (I l ) € IC0

+

+

(11.18)

In conclusion, the fictive charge surface density a-may be generated by multiplying ea+ multipolar component Gon(l,m) of the normal component of bo, (on the sphere of radius a ) by a factor AI) which only depends on I but not on m.

+ l)(€(I + l ) € + IC0

(21

f(l)=

€0)

(11.19)

3€, 2€,

-1

C,

(11.23)

€,

The reason for this apparent disagreement is that in fact Mertius et al. calculate 1 t , - 121.1COS8 1 €,- 1, no = - --= - -Go,( 1) (11.24) 47r cy r3 4 a €, '

so that the exact surface charge density = fMST(1)"O

u

is given by (11.25)

whereas we calculate u

= f(Z)BO,(l)

(11.26)

At the end, the two processes lead to the same value of

U.

Calculation of the Thermodynamic and Transport Properties of Aqueous Species at High Pressures and Temperatures: Dissoclation Constants for Supercrltical Alkali Metal Halldes at Temperatures from 400 to 800 "C and Pressures from 500 to 4000 bar Eric H. Oelkers and Harold C. Helgeson* Department of Geology and Geophysics, University of California, Berkeley, California 94720 (Received: June 15, 1987; In Final Form: October 8, 1987)

Dissociation constants from 400 to 800 O C at 500-4000 bar have been calculated for 14 aqueous alkali metal halides from experimental conductance measurements reported by Quist and Marshall' using the Shedlovsky equation and the law of mass action, together with independently calculated limiting equivalent conductances for the electrolytes generated from equations given by Oelkers and Helgeson? The requisite activity coefficients were calculated for neutral and charged species from the Setchtnow and extended Debye-Hiickel equations, respectively. Where direct comparison can be made, dissociation constants reported by Quist and Marshall3 and Dunn and Marshall4 fall within estimated uncertainties of corresponding values generated in the present study. The logarithms of the calculated dissociation constants range from - 4 . 7 at low temperatures and high pressures to --4.5 at high temperatures and low pressures.

Introduction Dissociation constants computed from experimental measurements of the conductances of supercritical aqueous electrolyte solutions have been reported by Fogo, Benson and Copland,' (1) Quist, A. S.; Marshall, W. L. J . Phys. Chem. 1969, 73, 987. (2) Oelkers, E. H.; Helgeson, H. C. Geochim. Cosmochim. Acra 1988,52, 63. (3) (a) Quist, A. S.; Marshall, W. L. J . Phys. Chem. 1968, 72, 684. (b) Quist, A. S . ; Marshall, W. L. J . Phys. Chem. 1968, 72, 2100.

0022-3654/88/2092-163 1$01.50/0

Franck? Quist et al.,' Quist and Marshall>* Ritzert and Franck: Dunn and M a r ~ h a l l Mangold ,~ and Franck,lo and Frantz and (4) Dum, L. A.; Marshall, W. L. J. Phys. Chem. 1969, 73, 723. (5) Fogo, J. K.; Benson, S . W.; Copeland,C. S . J. Chem. Phys. 1954,22, 212. (6) (a) Franck, E. U. Z . Phys. Chem. 1956,8,92. (b) Franck, E. U. Z . Phys. Chem. 1956,8, 107. (c) Franck, E. U. Z . Phys. Chem. 1956,8, 192. (7) Quist, A. S.; Franck, E. U.; Jolley, H. R.; Marshall, W. L. J . Phys. Chem. 1963, 67, 2453. (b) Quist, A. S.; Marshall, W. L.;Jolley, H. R. J. Phys. Chem. 1965, 69, 2726.

0 1988 American Chemical Society

Improvements of the continuum model. 1. Application to the calculation

Recommend Documents