J. Phys. Chem. 1994,98, 5147-5151
5147
On the Solubility of Aqueous Electrolytes John Perkyns and B. Montgomery Pettitt' Chemistry Department, University of Houston, Houston, Texas 77204-5641 Received: November 15, 1993; In Final Form: February 24, 1994"
The existence of a high-density phase separation for molecular and continuum solvent models of ion-water solutions are demonstrated using the DRISM/HNC theory. The molecular level effects which influence solubility and the inability of continuum models to exhibit some of these effects are discussed. The general trends in experimental alkali halide solubilities a r e found to be consistent with the present model calculations. The dependence of the phase separation on ion size is found to be primarily entropic in nature.
Introduction The last half-century has seen an enormous growth in the study of electrolyte solutions, yet by far the greatest proportion of effort has been expended on aqueous systems where the concentration ofsalt is below 1 M. While knowledge of the structureof relatively dilute solutions has been greatly enhanced by such attention, the structure and phenomena of higher concentrations is less well understood. Of long-standing interest is the difference in solubility between ions of the same charge magnitude.' The alkali halides are one of the best studied periodic examplesZ which show nontrivial trends.3 The patterns found in ionic solubilities are due to a combination of both the transfer free energies from the gas phase as well as lattice energies. The latter have been well studied, and quantitative models have long been known.4 More recently the gas-phase transfer part of the Born cycle of solvation has been considered by simulations for single ions at infinite dilution (see for instance ref 5). However, free energies of solvation at finite concentration and solubility limits have proven less tractable. Recent advances in integral equation theories of finite concentration aqueous ionic solutions have yielded good qualitative results as regards a variety of physical properties for both angular expansion6 and site-site model^.^ Using the distribution functions from these techniques in combination with a statistical thermodynamic theory of stability8 allows us to consider some of the characteristics of aqueous salt solutions near saturation using atom-based models. The difficulties in modeling what is occurring at a molecular level near the solubility limit of saturated salt solutions are numerous. First it is unlikely that a simple classical (meanfield) interaction potential model of a strong electrolyte will successfully duplicate experimental thermodynamic quantities for both the solvated and the solid salt. The most commonly used models for solvated ions have pair potentials which combine a Coulombic interaction with a Lennard-Jones 6-1 2 functional form9-'2 to model the short-ranged repulsive and dispersive forces. However, models of ionic crystals usually only share the Coulombic term, and often include as a short-ranged contribution an exponential repulsion such as in the Huggins-Mayerl3 form. Second, there is no molecular theory that describes both solvated and solid salt phases in equilibrium. Neither do there exist sufficiently accurate theories for either solvated or crystalline salt phases such that the thermodynamic quantities predicted can be relied on to duplicate the corresponding region of an experimental phase diagram. For computer experiments the time scale or sampling problems of the dissolution process provide nearly insurmountable difficulties at present. Nevertheless, a molecular description of ionic correlations and the solvation effects occurring near solubility limit conditions, @
Abstract published in Aduance ACS Abstracts, April 1, 1994.
even in the context of idealized models and approximate theory, could form a basis from which to describe a number of systems of current interest from phase equilibria and critical phenomena to the salting out of proteins. As yet the only theoretical treatment of phase separations in high concentration aqueous salt solutions are in the context of continuum theories for models such as the restricted primitive model (RPM).14 The phase stability phenomenon we will examine here is found at high density and temperature ( p , T ) on the system phase diagram. However, integral equations and simple ion models have been used to elucidate other phase-stability phenomena found elsewhere on the phase diagram. Since the existence of a lowdensity, low-temperature, ( p , T ) critical point in charged hardsphere fluids was shown,I5the onset of phase instability has been examined theoretically in several contexts. Analytical solutions of the hard-core Yukawa fluid (HCYF) in the mean-spherical approximation (MSA)lGISand for the adhesive hard-sphere fluid (AHSF) in the Percus-Yevick (PY) appro~imation'~J0have demonstrated a dependence on the combination of theory and model of important characteristics, such as whether the oftendemonstrated phase region where no numerical solution is possible coincides with the true spinodal line, or which values of critical exponents are found. For theory/model combinations where no analytical solution is possible, such as HNC/RPM or PY/RPM, the picture is not so clear, and the determination of the position of the true spinodal line, if it is exhibited, and the values of critical exponents, is extremely difficult.z1 The phase phenomenon we wish to elucidate is, however, quite distinct from that which is manifested in the low ( p , T ) (or equivalently low dielectric constant, t, and low density) region of the RPM. If solubility limits in aqueous electrolytes are to be understood, the high ( p , ~ region ) of the phase diagram must be examined. It has been shown that a phase separation in this region exists for a mixture of hard spheres with positively nonadditive diameters and that the phase separation as predicted by the H N C theory seems to demonstrate a true spinodal decomposition line coinciding with the line bounding the region of no numerical solution.z1 It is not possible to prove exact coincidence, but it may be concluded that the H N C theory in combination with this particular model is much closer to selfconsistency in this regard than for the low ( p , T ) phase separation. It is not our purpose here to demonstrate whether such spinodal lines are likely to be correctly predicted by integral equation theory but simply to show that for simple models of aqueous monovalent electrolytes, in both a continuum solvent/HNC treatment and in a full molecular solvent/HNC treatment, a high ( p , 7') region of thermodynamic instability exists and further that its dependence on ionic concentration is itself highly dependent on ion size. Nevertheless, the formal aspects of the theoretical predictions of phase stability for useful models are important to
0022-3654/94/2098-5 147%04.50/0 0 1994 American Chemical Society
The Journal of Physical Chemistry, Vol. 98, No. 19, 1994
5148
the future improvement of such theories, so we have examined stability conditions where appropriate.8 In this paper we present results for the 1-6-12 model of electrolytes in a dielectric continuum and a full molecular waterlike solvent. The molecular solvent theory used is the dielectrically consistent reference interaction site theory (DRISM) recently prop0sed.~,~2The maximum attainable salt concentration for which a solution to the theory can be found for a range of anion and cation sizes is presented, and the distinct features of the results are discussed. The vast differences between results for the two models are interpreted in terms of two distinct molecular effects: the tendency of ions to form pairs or aggregates, and a hydrophobic effect. While both effects are in competition for all ion sizes studied, the wide range of maximum attainable concentrations is dependent on which effect dominates.
Theory
Perkyns and Pettitt strictest criterion for the full molecular solvent salt solution model is that of compositional stability: ( a d d x i )~ (dPi/dXi)$
-
1
, p
VXJ
+ P,(G+- - 2G+u + GUJI
2 0 (7)
wherexiis themole fractionofspecies i,piis thechemical potential of species i, id denotes the quantity for a system behaving ideally, and Gij is the usual "Kirkwood G"25 defined by
G , = 4 s c r l h i j ( r i j )dr, Throughout we use the subscripts +, -, u, and u to denote cation, anion, solute, and solvent species, respectively, and for the salt MP+X-*-we have v = v+ v-. For the monovalent models used here Y = 2. For the related continuum model the strictest stability criterion is the mechanical criterion
+
The DRISM theory used here has been presented previously, and the reader is referred to ref 7 for details. All intermolecular pair potentials (between sites i and j ) are of the form
where uijand ci,are the usual diameter and well-depth parameters. For the molecular solvent the parameters for SPCE water were employed.23 The DRISM theory7J2 corrects the dielectric inconsistencies inherent in the RISM theory24by including "missing"corre1ations in the place of bridge functions, bij(r). We write
c,(r) = exp[-@u,(r)
+ tij(r)]- 1 - t,(r) + bij(r)
+
= cHiij(r) b,(r)
(2)
which defines cRij(r)as the solution to the usual H N C equation, and tij(r) = hij(r)- cHiij(r).Here hij(r)and cHiij(r)are the total and direct correlation functions, respectively, between sites i and j . Assuming all the functions bj(r) can be Fourier transformed, we can build matrices b, &,and bout of the Fourier space functions &ij(k), E ~ ; i j ( k )and , hij(k), respectively, and write the RISM equation in the form p&
= b(5H
+ b)b + b(& + b)p&
(3)
where is the matrix of intramolecular site correlations, and p is the matrix of site densities. Now defining z via ji = 580
+ 682
(4)
we can rewrite eq 3 as (p&
- x) = ( b + g)&(b
+ 2 ) + ( b + g)CH(p&?
- 2)
(5)
It has been shown7 that dielectric consistency can be achieved by equating the components of the matrix 2 to jiij(k>
= -jd-k4J jO(-kdiJ jl(-kdiz) X
where X T , ~is the isothermal compressibility at constant mole fraction. Before we discuss our results, some remarks about eq 7 are in order. The equality between the first and second terms in eq 7 is given in ref 26. It is derived from the extension of KirkwoodBuff theoryZ7 to solvent-electrolyte mixtures. We note that this equality is exact. That the partial derivativeof chemical potential with respect to mole fraction must be nonnegative, for components in an equilibrium mixture, is a well-known r e ~ u l t . ~ 8 -Whether 2~ it is the strictest condition for stability depends on the particular system under consideration. When Kirkwood G's from an approximate theory are applied in eq 7, a question arises. Are all points in the phase diagram complying with this criterion coincident with computationally available points? While exact agreement of these two methods of predicting the spinodal decomposition line is not proven, no discrepancy was discovered. The recent work of Ursenbach and Pateys which applies eq 7 to ions in a solvent with interactions modeled by angularly dependent potentials, discovers the same close agreement. However, the position on the phase diagram of the spinodal decomposition predicted by approximate theories is not exact for the respective models.
Results and Discussion Calculations for 1-1 electrolyte solutions with a range of values for u++ and u- have been performed using the present theory. For each calculation a solution was obtained for a low salt concentration at the desired combination of ion sizes, and then the concentration of salt was increased by small increments until either a 12.0 M solution was achieved or a concentration was found above which no solution would converge numerically, [M+X-I,,,. Nonphysical solutions were in some cases obtained but were clearly distinguished by large negative calculated dielectric constants and large deviations from the expected charge neutrality conditions:
jo(kdjx) jo(kdjy) jl(kdjZ)[ P & , ( ~ ) I (6a)
where j i ( x ) is the spherical Bessel function of order i and
JJk) =
["YPU
3 1 exp(-ak'),
a>0
(6b)
with pu the density of the solvent and y is the usual dipole density. DRISM theory results are found by solving eqs 5 and 2 simultaneously. The criteria used here to discover the concentration marking theonset of absolute instability for ionicsolutions were first applied to integral equation calculations by Ursenbach and Patey.8 The
There was no physically reasonable solution found where the stability criteria (eqs 7 or 9, where appropriate) were violated, and no false solution found where the stability criteria were reasonable. All false solutions diverged when the numerical range of integration of the correlation functions was doubled. All calculations were carried out at a temperature of 300 K and with constant ion-ion Lennard-Jones ti, values of 0.1 kcal/mol. Two unavoidable approximations which, when taken separately, become increasingly less satisfactory with increasing concentration
Solubility of Aqueous Electrolytes
The Journal of Physical Chemistry. Val. 98, No. 19. 1994 5149
5
3
%
4
%
ti
ti 2
3
2 1
2
a++(a)
3
4
Figure1. [M+X-], surface forthe molecular solvent model as a function of r.. and c++. The shaded areas above I2 0 M represent regions where valuesof [MIX-], wereconsidered tw high tobe reliable. Thedashed line outlines the shaded region where no solution could be achieved at any finite salt concentration
were used throughout, but when taken together they largely offset each other with regard to the concentration where a numerical solution ceases to be possible and do not affect the general trends in structure or interpretation. First, the concentration was increased by keeping the total number density constant (at that ofpure water.0.0333 molecule/A3), i.e.,each ion addedreplaced a water molecule. This approximates constant-pressureconditions and is reasonableonly for concentrations below about 1.0 M, but correct densities at I atm for the ion sizes necessary to examine the general physical processes involved with solvation areobviously not available. Second, the dielectric decrement associated with ionic strength as predicted by HNC-based theories7.26.30 also becomes more unsatisfactory with increased concentration. Sample calculations for NaCl were performed for which experimental (constant pressure) densities" and dielectric constants32 were available. As expected, the experimental dielectricconstants which decrease more rapidly than the H N C predictions, lowered the value of [Na+CI-I,, considerably when taken alone. The experimental densities, which represent a larger number of water molecules per ion than the constant total density approximation, increased thevalueof [Na+CI-],,considerably. However, when thetwoapproximations weretaken together, thevalue [Na+Cl-1, differed from the calculation using experimental input quantities by less than 10%when the total concentration was below 12 M. The values of [M*X-I... vs a range of cation and anion sizes for the molecular solvent are given in Figure I . The effect of the asymmetry of charge density around a water molecule manifests itself as an asymmetry with respect to ion size, particularly demonstrated by the positions where twoareas of high [M+X-I,. occur. By contrast, the continuum solvent results are simply a steeply increasing function of u+ (Figure 2). The complex structure of Figure 1 is due to two distinct molecular effects which will be examined in turn. Figure 1 contains a large region of very low solubility, corresponding to the case when both ion sizes are small, which results when cation-anion interactions overcome ionsolvent interactions, causing the formation of ion pairs or aggregates. This region is bounded by a narrow section of (u++,u..) values where [M+X-]...risesdramatically . Anexplanationis found by comparingthechangein thewelldepthoftheinterionicpotential, uc(r) (Figure3),and thechangein thedepthofthe first minimum in the interionic potential of mean force (PMF), o+.(r)(Figure 4). for the change from a (c++ = 1.0 A, u.. = 4.0 A) solution to a (u++ = 1.0 A, u.. = 4.5 A) solution. The relatively small %decrease in potential well depth, which occurs when the anion is increased in size by 12%. has resulted in a 35% decrease in the first minimum of the corresponding PMF (while keeping the
1 2
1
a++(a)
3
FIgureZ. [M+X-],surfacefor thecontimumsolvent modelasafunction of v.. and st+. The shaded areas are as in Figure I .
0
2
r(i)
4
6
Figure 3. Pair potentials for oppositely charged ions (0 = I fksT, ks the
Boltimamconstant). Thesolidand dashed lines represent the molecular model tentials for ions of sizes (o++ = 1.0 A, e. = 4.0 A) and (c++ = 1.0 fo.. = 4.5 A), respectively. The small circle line represents the continuum solvent potential for the latter, with a continuum e = 80.
1
r(A)
5
7
Figure 4. Potentials of mean force for oppositely charged ions. The solid, dashed, and dot-dashed lines represent the systems (v++ = 1.0 A. 0.. = 4.0 A, [M+X-] = 0.095 M), (c++ = 1.0 A. m. = 4.5 A, [M+X-] = 0.095 M) and (u++ = 1.0 A,..u = 4.5 A, [M*X-] = 20.0 M).
respectively. The small circle line corresponds to the continuum model for the system ($++ = 1.0 A, c.. = 4.0 A, [MIX-] = 0.095 M). ionic concentration constant). Thus a fine balance between the strength of interionic f o r m and ion-solvent formcauses dramatic differences between the degree of solvation of ionic solutes. The effect of change in anion ion size on the distribution of cations is shown in Figure 5 , where a huge peak in g++(r) representing correlations between positive ions indicates the presence of ion pairs or aggregates. The peak position indicates that pairs ofcationsare tooclose tobesolvent-or anion-separated
5150
The Journal of Physical Chemistry, Vol. 98, No. 19, 1994
6
Perkyns and Pettitt
01
I
I
I
1
4.5
5.0
4 g++w
2
0
1
3 r(A)
5
3.0
7
Figures. Radialdistributionfunctionsfor cation pairs. Thelinesrepresent the same systems as in Figure 4.
pairs, while not close enough to be in direct contact. The restructuring of the solution from high to low probability of finding two cations in proximity is consistent with the separation of ion pairs (Figure 4). The dramatic drop in this peak to about 15% of its original value when the opposite (anion) ion size is changed from u-- = 4.0 8, to u-- = 4.5 8, is a purely indirect effect of changing the counterion size. The pair potential between cations and the ionicconcentration are unchanged for these two functions. The significant difference in structure between these two solutions is accompanied by a much smaller than expected change in average total excess internal energy per particle, B( U )/N, which changes (becomes less favorable) by only 0.1%. A componentwise breakdown of contributions to the internal energy reveals almost exactly offsetting changes in the ion-solvent and ion-ion terms, while the solvent-solvent contribution changes by less than 0.01%. This system restructuring must therefore be accompanied with a state of more favorable entropy. The effect of increasing the concentration of salt from0.095 to 20 M is also included in Figures 4 and 5, to demonstrate that the effect of increasing the concentration by a factor of over 200 has less effect on the overall structure of the system than changing the size of one ion by 12% in this particular region of solution space. This effect of small ions pairing and forming aggregates is also manifested in the values of [M+X-Imax as shown in Figure 2 for the continuum model. The ion-ion correlations given by this model are also given in Figures 4 and 5 , but the relative shallowness of the corresponding pair potential (Figure 3) gives rise to fairly structureless distributions. The fine balance between the dominance of solvated and paired ions is not present in such a simple potential model, so the decrease in [M+X-],,, is not as steep, but the general trend of lower [M+X-I,,, values for smaller ions is still observed. Density contours for this model are in general shifted to considerably smaller ion sizes than those for the molecular solvent model. This is primarily because this model potential is orders of magnitude too weak at short range (Figure 3) with the corresponding lack of solvation shell structure induced in its interionic correlations (Figures 4 and 5 ) . Figure 1 also shows a decrease in [M+X-],,, as either ion diameter is increased. To elucidate this effect, a breakdown into components of the contributions to-@( U)/Nfrom the interactions between different species types is given in Figures 6 and 7, isolating the increase in just the anion size while the cation size and total salt concentration are held constant. Of most interest here is the decrease in the anion-solvent contribution to the internal energy, coupled with the simultaneous increase in the ion-ion and solventsolvent contributions. The increase in configurational energy among like species with a concurrent decrease between unlike species is consistent with the hydrophobic effect. This can be understood in terms of the changing strengths of the electrostatic potentials as ion sizes increase. The strength of the electric field
3.5
4.0
.--(A) Figure 6. Contributions to the average excess configurationalenergy per particle in the system for anion sizes ranging from u-- = 3.0 A to u-= 5.0 A. The solid, dot-dashed, and dashed lines represent cation-
solvent, ion-ion, and anion-solvent contributions, respectively. The salt concentration and cation size are fixed at 0.25 M and u++ = 4.0 A, respectively.
-38
/? N -39
-40 3.0
3.5
4.0
4.5
5.0
.--(A) Figure 7. Total (small circle line) and solvent-solvent contribution (dot-
dashed line) to the average excess configurational energy per particle. The values are for the same systems as Figure 6 . at the “surface” of a (fixed charge) ion becomes lower as the size of the ion increases. The ability for an ion to restructure the solvent in its vicinity is therefore decreased, and the structure of the water near the surface of the ion becomes more like that of the bulk, with a corresponding increase (to a more negative value) in the average solvent-solvent energy. A molecular property sometimes associated with the hydrophobic effect is that of dewetting. The position of the maxima in the first solvation shell peaks in the gvi(r) functions would be expected to move to a larger r value than can be accounted for by the increase in the size of the ions alone if dewetting were present, but such an effect is not observed in these calculations. The first shell maximum in g,,,-,(r) for the water oxygen sites might also be expected to move to a slightly lower r value, but this is not seen either. The differences in the contributions to the internal energy must therefore be primarily due to rotational rearrangement of water molecules. This is entirely consistent with the almost constant total system configurational energy for the calculations represented in Figures 6 and 7. Since the contribution to the system free energy from the internal configurational energy is almost constant as the size of the ions is changed, any variations in the free energy, and hence the phase stability of the system, must be primarily entropic in nature. The two regions of highest [M+X-Imaxoccur (Figure 1) when one ion is relatively large and the other relatively small, with a saddle point at approximately (a++= 2.7 A, u-- = 4.3 A) between them. When one ion is small enough to form a tightly held solvation shell, yet the other is large enough that its interionic
Solubility of Aqueous Electrolytes attraction is not strong enough to penetrate into the shell, regions of strong solvation result since neither of the effects discussed above which favor phase separation can occur. The saddle point between these high [M+X-Imaxregions results when increasing thesmaller ionsizedecreases the strength with which thesolvation shell is held, allowing the onset of the restructuring of water molecules around it leading to gradual competition by the hydrophobic effect as discussed above. Unfortunately, comparison of [M+X-I,,, to experimental solubilities (even if values for all monovalent salts were available at the same thermodynamic state) is inappropriate since the inclusion of other necessary factors such as the free energy of formation of the appropriate salt crystals are beyond the scope of these calculations. However, general trends in alkali halide solubilities3 are consistent with the general trends discovered here. When both ion sizes are small as in the case of LiF, the solubility in aqueous solution is also low, 0.10 M a t 18 OC. When both ions are increased to relatively large sizes, as in the change from RbBr to CsI, the solubility limit also decreases (5.93 M at 5 OC and 1.69 M at 0 "C). And last, when one ion is large and the other small, as in CsF and LiI, the solubility limits are relatively high (25.16 M a t 18 OC and 12.4 M a t 20 OC, respectively). Furthermore, thereseen to be no solubility limits among the alkali halides which break these general trends. Summary and Conclusions The existence of a high concentration phase separation in both a full molecular solvent model of a salt/water mixture and a continuum solvent model has been demonstrated for our site-site HNC-based integral equation theory. The maximum salt concentrations attainable using this theory were calculated as a function of ion size. It was discovered that in the case of the molecular solvent, the values of [M+X-],,, are decreased by each of two molecular effects, the tendency of small ions to form pairs and aggregates, and the onset of the hydrophobic effect as both ions are increased in size. The corresponding continuum model only displayed the former effect, at considerably lower values of ion size. Changes in the various contributions to the average excess configurational internal energy were found to almost exactly offset each other even when the molecular distributions were substantially altered by changes in ion size, leading to the conclusion that the structure of the [M+X-Imaxsurface is largely entropic in nature. It is well-known that unlike ion pairing is also entropically dominated.33 This entropic effect is due to the temperature derivative of the dielectric constant of water. The general trends exhibited by such a [M+X-I,,, surface for this solvent model and theory are entirely consistent with the trends in experimental solubility limits.
The Journal of Physical Chemistry, Vol. 98, No. 19, 1994 5151 Acknowledgment. The authors thank the Robert A. Welch Foundation, the N.I.H., the N.S.F., and the Texas Coordinating Board for partial support of this work. The N.S.F. is thanked for local computational support. B.M.P. acknowledges the Alfred P. Sloan foundation for support. The authors thankc. Ursonbach and Prof. G. Patey for a preprint of ref 8. References and Notes (!) Conway, B. E. Ionic Hydration in Chemistry and Biophysics; Elsevier: New York, 1981. (2) Burgess, J. Ions in Solution: Basicprinciples of chemical interactions; Ellis Harwood Ltd.: New York, 1988. (3) Weast, R. C., Ed. CRC Handbood of Chemistry and Physics, 66th ed.; CRC Press, Inc.: Boca Raton, FL, 1986. (4) Tosi, M. P.; Fumi, F. G. J. Phys. Chem. Solids 1964, 25, 45. (5) Straatsma, T.; McCammon, J. A. Annu. Reu. Phys. Chem. 1992,43, 407. (6) Kusalik, P. G.; Patey, G. N. J . Chem. Phys. 1988, 88, 7715. (7) Perkyns, J. S.; Pettitt, B. M. J . Chem. Phys. 1992, 97, 7656. (8) Ursenbach, C.; Patey, G. N. J . Chem. Phys. 1994, 100, 3827. (9) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swammathan, S.;Karplus, M. J . Comput. Chem. 1983, 4, 187. (10) Weiner, S.J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J . Comput. Chem. 1986, 7 , 230. (11) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. SOC.1988, 110, 1657. (12) Pranata, J.; Wierschke, S.G.; Jorgensen, W. L. J . Am. Chem. Sot. 1991, 113, 2810. (13) Huggins, M. L.; Mayer, J. E. J . Chem. Phys. 1933, 1 , 643. (14) Caccamo, C.; Malescio, G. J . Chem. Phys. 1989, 90, 1091. (15) Stell, G.; Wu, K. C.; Larsen, B. Phys. Rev. Lett. 1976, 37, 1369. (16) Cummings, P. T.; Smith, E. R. Chem. Phys. 1979, 42, 241. (17) Cummings, P. T.; Stell, G. J. Chem. Phys. 1983, 78, 1917. (18) Cummings, P. T.; Monson, P. A. J . Chem. Phys. 1985, 82, 4303. (19) Baxter, R. J. J. Chem. Phys. 1968, 49, 2770. (20) Fishman, S.;Fisher, M. E. Physica A 1981, 108, 1. (21) Belloni, L. J . Chem. Phys. 1993, 98, 8080. (22) Perkvns, J. S.;Pettitt, B. M. Chem. Phys. Lett. 1992, 190, 626. (23) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J . Phys. Chem. 1987, 91, 6269. (24) Chandler, D.; Andersen, H. C. J . Chem. Phys. 1972, 57, 1930. (25) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic Press Inc.: London, 1986. (26) Kusalik, P. G.; Patey, G. N. J . Chem. Phys. 1987, 86, 5110. (27) Kirkwood, J. G.; Buff, F. P. J . Chem. Phys. 1951, 19, 774. (28) Prigogine, I.; Defay, R. Chemical Thermodynamics;LongmansGreen and Co.: London, 1954 (translated by D. H.Everett). (29) Tisza, L. Generalized Thermodynamics; MIT Press: Cambridge, MA, 1966. (30) Rasaiah, J. C.; Isbister, D. J.; Stell, G. J . Chem. Phys. 1981, 75, 4707. (31) Washburn, E. W., Ed. International Critical Tables; McGraw-Hill: New York, 1926. (32) Behret, H.; Schmithals, F.; Barthel, J. 2.Phys. Chem. (Munich) 1975, 96, 73. (33) Pettitt, B. M.; Rossky, P. J. J . Chem. Phys. 1986, 84, 5836.
