4218
J. Phys. Chem. B 2000, 104, 4218-4221
Macroion-Macroion Correlations in Presence of Divalent Counterions B. Hribar and V. Vlachy* Faculty of Chemistry and Chemical Technology, UniVersity of Ljubljana, Asˇkercˇ eVa 5, 1000 Ljubljana, SloVenia ReceiVed: December 13, 1999
Monte Carlo results for a primitive model of a micellar solution are presented. The macroions are assumed to carry 20 negative charges and the counterions are monovalent and/or divalent. The ratio of diameters of the ionic species, described as charged spheres, is 3.0 nm:0.4 nm (macroions:small ions). An important conclusion of the study is that by replacement of monovalent counterions with divalent ones it is possible to change the structure of the solution. In solutions with predominantly divalent counterions, the macroions approach each other and share the layer of counterions. This structure is rather stable; in our case the macroion clusters do not “dissolve” upon dilution from 0.02 to 0.01 mole dm-3 and are not significantly affected by moderate addition of simple electrolyte. Thermodynamic considerations indicate an energy-entropy compensation in the process of dilution.
1. Introduction Aqueous solutions of charged colloids, surfactant micelles, and polyelectrolytes play an important role in our everyday life. Detergents, dyes, and globular proteins are only a few obvious examples. Because of their interesting properties and widespread application, these systems represent a challenge for both experimentalists and theoretical scientists. Experimental evidence suggests that charged macroions are stabilized by electrostatic forces, so it is assumed that the DLVO theory1 is applicable to these systems. This classical theory assumes competition between the electrostatic repulsion of an overlapping electrical double layer and the short-range van der Waals attraction. In recent years experimental results have been reported which are not consistent with the classical DLVO approach.2-6 Reports about the occurrence of large stable voids in otherwise homogeneous suspensions, equilibrium separation into two phases corresponding to a liquid-vapor phase transition, and metastable colloidal crystallites (the latter results apply to confined geometry)4,5 suggest the existence of a strong attractive interaction of Coulombic origin acting between similarly charged colloidal particles. Note that this cannot be the van der Waals type of attraction,1 for the latter is short ranged. Can like-charged macroions in solution attract each other? The idea of a net attraction of electrostatic origin, which can happen due to the presence of small ions (counterions and possibly co-ions) in the system, has been voiced in the scientific community for quite some time7,8 (for reviews see refs 9 and 10). Despite different explanations how such an attractive force may arise,11-16 there is concern that some of the results are simply artifacts of the measurements.17 The most compelling evidence that the interaction between equally charged surfaces can be attractive, comes from measurements of the force between two charged mica surfaces immersed in an electrolyte with divalent counterions.18,19 Moreover, these experimental results are supported by theoretical predictions. Of course, an * Corresponding author. Phone/Fax: +386-61-125-44-58. E-mail:
[email protected].
effective interaction between mobile spherical macroions is different than that for two infinite charged planes, and this problem requires a separate investigation. In a previous paper20 we published Monte Carlo results for -20:+1 and -20:+2 asymmetric electrolytes in the primitive model approximation. The results for the macroion-macroion pair distribution function (pdf), presented in ref 20, indicate the important role of correlations between the counterions in solution. For a model solution with monovalent counterions, the position of the peak of the macroion-macroion pdf is shifted toward larger distances on dilution, as expected from the repulsive interaction between the bare macroions. On the other hand, for a solution with divalent counterions (-20:+2 electrolyte), the peak of the pdf does not change its position upon dilution. The two macroions stay close to each other, sharing a layer of counterions. This counter-intuitive result (the ions in ref 20 are modeled as charged hard spheres with no short-range attractive force acting between them) prompted a more systematic study of the role of interionic correlations in solution, presented in this paper. In particular, we were interested to see how the structure of the solution gradually changes by replacing monovalent counterions with divalent ones. The system under study is therefore a solution of macroions with a mixture of mono- and divalent counterions. 2. Model and Methods The ions were treated as charged hard spheres and the solvent was considered as a continuum with a dielectric constant r equal to that of bulk water at T ) 298 K. In this model, the interaction potential between two particles of species a and b is given by
{
r < (ra + rb) ∞, 2 uab(r) ) e zazb r g (ra + rb) 4π0rr,
(1)
Here, r is the distance between the particles and ra is the radius of a particle of type a. Also, za is the valency of an ion of type a and e is the proton charge. The indices a and b stand for
10.1021/jp9943732 CCC: $19.00 © 2000 American Chemical Society Published on Web 04/06/2000
Macroion-Macroion Correlations
J. Phys. Chem. B, Vol. 104, No. 17, 2000 4219
macroions (p), counterions (c), and co-ions (k). The radii of the ions were rp ) 1.5 nm and rc ) rk ) 0.2 nm, and the corresponding valencies were zp ) -20, zc ) 1, or zc ) 2 and zk ) -1. Most of the results presented here were obtained by the Monte Carlo method;20-23 in some cases the simulation data were supplemented by the hypernetted chain integral equation calculations (HNC).22,23 The reason for not using the HNC approach in all cases lies in fact that this theory yields convergent results only for a limited range of model parameters.23 The HNC approximation consists of two equations. One is the Ornstein-Zernike equation, which correlates the total correlation function, hij(12) ) gij(12) -1, and the direct correlation function, cij(12)24
hij(12) ) cij(12) +
∑k Fk ∫cik(13)hkj(32)d(3)
(2)
where Fk is the number density of the component k. The second equation is given by 24
cij(12) ) hij(12) - [1 + hij(12)] - βuij(12) + Bij(12)
(3)
where β is 1/kBT; T is the absolute temperature and kB Boltzmann’s constant. The HNC approximation sets the unknown bridge function Bij(r) to zero for all distances r. A necessary first step toward numerical solution of the OrnsteinZernike equation for ionic systems is the renormalization procedure. For the three-dimensional ionic fluids, the renormalization scheme is well known (see for example, 23) and will not be repeated here. Computer simulations do not contain the usual statisticalmechanical approximations, and they provide exact results for a given physical model.25 In this way the results presented here are free of the assumptions of the theories11-14,16 designed to explain deviations from the DLVO theory. The Monte Carlo simulations were performed in a canonical ensemble (at constant volume and temperature) with 64 (or in some cases 128) macroions (p) and an equivalent number of counterions (c) and co-ions (k) in the system. Long runs are needed to obtain reliable results: the averages were taken over 50 to 70 million configurations with at least 5 million configurations spent for the equilibration. To avoid effects due to the finite size of the simulation cell, the Ewald summation method25 was used. The excess internal energy per particle was calculated by
Eex N
)
2π Ft
FiFj ∫∞0 uij(r) gij(r)r2 dr ∑ i,j
(4)
where Ft is total number density of particles. The osmotic coefficient, defined as the ratio of the calculated pressure and ideal pressure Φ ) P/Pid, was obtained via the virial route
Φ)1-
2 πβ 3 Ft
FiFj ∫∞0 u'ij(r) gij(r)r3 dr + ∑ i,j 2π
FiFj(ri + rj)3 gij(rij) ∑ 3F i,j
(5)
t
where gij(rij) is the value of the pair distribution function in contact, rij ) ri + rj. More numerical details concerning the Monte Carlo method and HNC calculations are given in references 21-23. 3. Results and Discussion Distribution Functions. We begin this discussion by considering the macroion-macroion pair distribution functions,
Figure 1. Pair distribution functions for macroions at different values of x (cp ) 0.02 mole dm-3) as obtained by the Monte Carlo simulation. The continuous curves in the far right and far left are obtained for pure -20:+1 (x ) 0) and -20:+2 solutions (x ) 1.0), respectively.19 The symbols denote results for x ) 0.1 and the dashed line denotes the pdf for x ) 0.5.
gpp(r), which are presented in Figure 1. An important variable of the present study is the degree of neutralization of the macroions with divalent counterions, x, defined as
x)
2c2 (c1 + 2c2)
(6)
where c2 and c1 are the concentrations of divalent and monovalent counterions, respectively. For x ) 0 (x ) 1) we have only monovalent (divalent) counterions present in solution, and x ) 0.5 means that one-half of the charge on the macroions is neutralized by divalent counterions. The results presented here apply to cp ) 0.02 mole dm-3. The height of the peak, which reflects ordering in the solution, does not change in a monotonic fashion when the monovalent counterions are replaced by divalent ones, but rather it passes through a minimum. As we can see from Figure 1, the height of the peak of the macroionmacroion pair distribution function at x ) 0.5 (dashed line) is smaller than the corresponding values at x ) 0 and at x ) 1. The results for x ) 0 (continuous curve on the right) and x ) 1 (continuous curve on the left), i.e., for pure -20:+1 and -20: +2 solutions, have been discussed before20 and are here presented for the sake of completeness. The pair distribution function denoted by symbols applies to x ) 0.1. By increasing the fraction of divalent counterions in the solution, the position of the first peak in the macroion-macroion pdf is changed; the macroions begin to approach each other. This is shown in Figure 2, where the position of the peak of this pdf is presented as a function of x at cp ) 0.02 mole dm-3. As we can see, the position decreases monotonically from 4.42 nm at x ) 0, until about 95% of the charge is neutralized by divalent counterions and remains constant after that. The broken line represents the HNC results, and symbols are Monte Carlo data. An influence of addition of simple electrolyte to polyelectrolyte solution is of substantial interest. In Figure 3 we show the macroion-macroion pair distribution functions in the absence (symbols) and presence (line) of added simple electrolyte. These results apply to -20:+2 electrolyte at cp ) 0.02 mole dm-3, while the concentration of added -1:+2 electrolyte is ck ) 0.02 mole dm-3. Note, that ck is the concentration of
4220 J. Phys. Chem. B, Vol. 104, No. 17, 2000
Hribar and Vlachy
Figure 2. Position of the first peak of the macroion-macroion pdf as a function of x (cp ) 0.02 mole dm-3). Symbols denote Monte Carlo data and the dashed line represents the HNC results.
Figure 4. The concentration dependence of the osmotic coefficient at three different values of x. The triangles apply to a -20:+1 solution (upper curve, x ) 0) and the circles to a -20:+2 solution (x ) 1.0). The symbols (squares) in the middle are Monte Carlo data for x ) 0.5. The results represented by lines were obtained by the HNC integral equation.
TABLE 2: Internal Energy E/NkBT, Entropy S/NkB, and Free Energy A/NkBT, Changes upon Dilution from cp) 0.02 mole dm-3 to 0.01 mole dm-3
Figure 3. Pair distribution functions for macroions in pure -20:+2 solution (line) and in mixture with added -1:+2 electrolyte (symbols); cp ) 0.02 mole dm-3 and ck ) 0.02 mole dm-3. The results are obtained by the Monte Carlo simulation.
TABLE 1: Excess Internal Energy, -Eex/NkBT, and Osmotic Coefficient, Φ, for Mixtures of -20:+1 and -20:+2 Electrolytes at Macroion Concentration cp ) 0.02 mole dm-3 x
-Eex/NkBT
Φ
0.0 0.1 0.3 0.5 0.6 0.75 0.8 0.9 0.95 1.0
3.929 ( 0.001 4.20 ( 0.01 4.837 ( 0.002 5.631 ( 0.002 6.108 ( 0.004 6.934 ( 0.004 7.258 ( 0.004 8.003 ( 0.005 8.382 ( 0.003 8.862 ( 0.004
0.85 ( 0.01 0.81 ( 0.01 0.75 ( 0.01 0.68 ( 0.01 0.64 ( 0.02 0.55 ( 0.01 0.52 ( 0.02 0.43 ( 0.01 0.43 ( 0.01 0.38 ( 0.02
co-ions in the mixture. As we can see from this figure, the two distribution functions almost coincide. This indicates that ionion correlations, which seem to be responsible for the net “attraction” between the macroions, are not affected by the presence of simple electrolyte. Thermodynamic Properties. The thermodynamic properties of these solutions may be of interest to many readers. In Table 1 the excess internal energy and osmotic coefficient Φ are presented as a function of x at concentration cp ) 0.02 mole dm-3. The concentration dependence of the osmotic coefficient
x
∆E/NkBT
∆S/NkB
∆A/NkBT
0 0.5 1
0.301 ( 0.002 0.313 ( 0.005 0.37 ( 0.01
0.627 ( 0.004 0.55 ( 0.01 0.40 ( 0.02
-0.326 ( 0.002 -0.233 ( 0.005 -0.02 ( 0.01
is shown in Figure 4 for three different values of x. The middle curve is calculated for x ) 0.5. The low value of the osmotic coefficient is mostly due to the strong binding of counterions in the field of macroions; for solutions with divalent counterions it is additionally decreased by “clustering” of macroions. The HNC results, wherever a numerical solution was possible, are also shown in this figure. In our previous study (cf Figure 2 of ref 20) we noticed that the peak of the macroion-macroion pdf in the -20:+2 solution does not change its position on dilution from cp ) 0.02 to 0.01 mole dm-3. We wished to understand this result in terms of the energy-entropy competition, reflected in the free energy difference between the final (diluted) and initial (concentrated) state of solution. Here we report data for the free energy (A), energy (E), and entropy (S) differences when the model solution is diluted from cp ) 0.02 to 0.01 mole dm-3. The free energy was obtained by integration of the Gibbs-Helmholtz equation in the form
β0A ) β0Aid +
∫ββ ) 0 Eex(β) dβ 0
(7)
where Aid represents an ideal part of the free energy and β0 applies to T ) 298 K. For the -20:+1 solution we obtained the result -406 J mole-1, while for the -20:+2 system we got a much smaller value of the free energy difference equal to -26 J mole -1. The results for ∆E, T∆S and ∆A are for x ) 0, x ) 0.5, and x ) 1.0 shown in Table 2. The thermodynamic properties of polyelectrolyte (-20:+2) electrolyte (-1:+2) mixtures are collected in Table 3 as a function of concentration of co-ions ck. We can see from these data that an addition of the -1:+2 electrolyte raises the osmotic coefficient and increases the excess internal energy of the solution.
Macroion-Macroion Correlations
J. Phys. Chem. B, Vol. 104, No. 17, 2000 4221
TABLE 3: Excess Internal Energy, -Eex/NkBT, and Osmotic Coefficient, Φ, for Electrolyte at cp ) 0.02 mole dm-3 with Added -1:+2 Simple Electrolyte of Different Concentration (ck is concentration of co-ions) -Eex/NkBT
Φ
ck [M]
-20:+1
-20:+2
-20:+1
-20:+2
0.0 0.0025 0.005 0.02
3.929( 0.001 3.891( 0.002 3.854( 0.001 3.647( 0.001
8.862 ( 0.004 8.697 ( 0.008 8.572 ( 0.004 7.910( 0.008
0.85 ( 0.01 0.85 ( 0.01 0.85 ( 0.01 0.88( 0.01
0.38 ( 0.02 0.37 ( 0.02 0.37 ( 0.05 0.43( 0.03
diseases, is connected with the separation of biological solutions into coexisting protein-rich and protein-poor phases.29 Acknowledgment. This work was supported by the Slovene Ministry of Science and Technology. Note Added in Proof. While this manuscript was under review, an important paper of V. Lobaskin and P. Linse30 was published. Their result about clustering of macroions in solutions with trivalent counterions is consistent with our findings (see also ref 26).
4. Conclusions
References and Notes
In this paper we complement our previous results20 and present further evidence of the importance of the correlation between divalent counterions in polyelectrolyte solutions. Our research indicates that in solutions with monovalent counterions the macroions are distributed at larger distances, while in solutions with predominantly divalent counterions they come closer to each other. The strong correlation between multivalent counterions (for example, if trivalent counterions are present) may even yield clustering of macroions and a nonuniform distribution with large voids in the solution.26 Our report20 about the existence of an attractive force between like-charge macroions was recently confirmed by Wu and co-workers,27 who studied the force between two isolated spherical macroions (infinitely dilute solution) in the presence of divalent counterions. Thermodynamic considerations indicate that for solutions with divalent counterions ∆E and T∆S contributions to ∆A of dilution from cp ) 0.02 to 0.01 mole dm-3 roughly cancel each other, yielding a relative small value of the free energy change ∆A. The computer simulations presented in this study are based on a well-defined model and are free of the assumptions of the ad hoc theories designed to explain disagreements with the DLVO approach. The conclusion is that the effective potential between the macroions depends on the composition and concentration of the solution in a more complicated way than suggested by the DLVO theory.1 The reason for the failure of the latter approach is, as discussed in more detail elsewhere9,15, the use of the Poisson-Boltzmann equation to describe the potential distribution around an isolated macroion, and on that basis, use of the superposition approximation for the macroionmacroion interaction. Our findings may have some implications for science and technology as well (see, for example, discussion in ref 28). In particular, knowledge about the stability of colloidal suspensions and polyelectrolyte solutions is crucial for industrial applications, and also for the biological sciences; a wide class of human diseases, known as protein condensation
(1) Verwey, E. J. W.; Overbeek, J. Th. G. Theory of the Stability of Lyophobic Colloids; Elsevier: New York, 1948. (2) Ito, K.; Yoshida, H.; Ise, N. Science 1994, 263, 66. (3) Tata, B. V. R.; Rajalakshmi, M.; Arora, A. K. Phys. ReV. Lett. 1992, 69, 3778. (4) Larsen, A. E.; Grier, D. G. Nature 1997, 385, 230. (5) Wickman, H. H.; Korley, J. N. Nature 1998, 393, 445. (6) Matsuoka, H.; Harada, T.; Kago, K.; Yamaoka, H. Langmuir 1996, 12, 5588. (7) Kirkwood, J. G.; Shumaker, J. B. Proc. Natl. Acad. Sci. U.S.A. 1952, 38, 863. (8) Oosawa, F. Biopolymers 1968, 6, 1633. (9) Vlachy, V. Annu. ReV. Phys. Chem. 1999, 50, 145. (10) Schmitz, K. S. Macroions in Solution and Colloidal Suspension; VCH: New York, 1993. (11) Ray, J.; Manning, G. S. Langmuir 1994, 10, 2450. (12) Lozada-Cassou, M.; Olivares, W.; Sulbaran, B.; Jiang, Y. Physica A 1996, 231, 197. (13) Sogami, I. S. Phys. Lett. A 1983, 96, 199. (14) Sogami, I. S.; Shinohara, T.; Smalley, M. V. Mol. Phys. 1992, 76, 1. (15) Kjellander, R. Ber. Bunsen-Ges. Phys. Chem. 1996, 100, 894. (16) Bowen, W. R.; Shariff, A. O. Nature 1998, 393, 663. (17) Palberg, T.; Wu¨rth, M. Phys. ReV. Lett. 1994, 72, 786. (18) Kjellander, R.; Marcelja, S.; Pashley, R. M.; Quirk, J. P. J. Chem. Phys. 1990, 92, 4399. (19) Kekichef, P.; Marcelja, S.; Senden, T. J.; Shubin, V. E. J. Chem. Phys. 1993, 99, 6098. (20) Hribar, B.; Vlachy, V. J. Phys. Chem. B 1997, 101, 3457. (21) Rebolj, N.; Kristl, J.; Kalyuzhnyi, Yu. V.; Vlachy, V. Langmuir 1997, 13, 3646. (22) Hribar, B.; Krienke, H.; Kalyuzhnyi, Yu. V.; Vlachy, V. J. Mol. Liq. 1997, 73, 74, 277. (23) Vlachy, V.; Marshall, C. H.; Haymet, A. D. J. J. Am. Chem. Soc. 1989, 111, 4160. (24) Friedman, H. L. A Course in Statistical Mechanics; Prentice Hall: Englewood Cliffs, 1985. (25) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University: New York, 1989. (26) Hribar, B.; Vlachy, V. Biophys. J 2000, 78, 694. (27) Wu, J.; Bratko, D.; Prausnitz, J. M. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 15169. (28) Murray, C. A. Nature 1997, 385, 203. (29) Liu, C.; Lomakin, A.; Thurston, G. M.; Hayden, D.; Pande, A.; Pande, J.; Ogun, O.; Asherie, N.; Benedek, G. B. J. Phys. Chem. 1995, 99, 454. (30) Lobaskin, V.; Linse, P. Phys. ReV. Lett. 1999, 83, 4208.
