Langmuir 1998, 14, 5077-5082
5077
Correlations between the Stability of Charged Interfaces and Ionic Exchange Capacity: A Monte Carlo Study Alfred Delville,*,† Najmeddine Gasmi,† Roland J. M. Pellenq,† Jean Michel Caillol,‡ and Henri Van Damme† Centre de Recherche sur la Matie` re Divise´ e, CNRS, 1B rue de la Fe´ rollerie, 45071 Orle´ ans Cedex 02, France, and MAPMO, De´ partement de Mathe´ matiques, Universite´ d’Orle´ ans, 45067 Orle´ ans Cedex 02, France Received March 10, 1998. In Final Form: June 24, 1998 (N,V,T) Monte Carlo simulations are performed by using hyperspherical geometry in order to investigate the mechanical stability and the ionic affinity of clay particles neutralized by a mixture of counterions of different valency. In this study, clay particles are described as uniformly charged surfaces and the interactions between ions and these charged surfaces are described within the context of the ‘primitive model’. The purpose of this study is to derive the maximum amount of monovalent counterions necessary to remove the stability of clay materials neutralized by divalent or trivalent counterions. From the equation of state of these charged materials we quantify the influence of ionic exchange on the interfacial free energy.
I. Introduction Suspensions of clay particles are the subject of numerous experimental and theoretical studies because of their use in many industrial applications (nuclear waste management,1,2 drilling,3 heterogeneous catalysis4). Their physicochemical properties (large adsorbing power, ionic exchange capacity, mechanical stability, acidity)3-6 result from the strong electrostatic couplings between the clay surfaces and their ionic, polar, or polarizable intercalates. Moreover, the understanding of the stability of clay materials at a fundamental level is crucial for various applications, such as the drilling of shales7 or the storing of ionic nuclear waste inside the interlamellar space of clays.1,2 Numerous numerical and theoretical studies have been devoted to this problem. The electric charges of natural and synthetic clays are generally neutralized by a mixture of counterions, while theoretical treatments used to analyze and predict the behavior of these charged interfaces generally exploit idealized models describing only homoionic materials.8-10 This oversimplification is not justified since the mechanical stability (cohesion or swelling) of these charged materials varies as a function of the valency of their neutralizing counterions.10-13 For instance, sodium montmorillonite is the prototype for an † ‡
Centre de Recherche sur la Matie`re Divise´e, CNRS. Universite´ d’Orle´ans.
(1) Adeleye, S. A.; Rautiu, R.; White, D. A.; Clay, P. G. J. Mater. Sci. 1995, 30, 583-586. (2) Olguin, M. T.; Solache, M.; Iturbe, J. L.; Bosch, P.; Bulbulian, S. Sep. Sci. Technol. 1996, 31, 2021-2044. (3) Manning, D. A. C.; Hall, P. L.; Hughes, C. R. Geochemistry of Clay-Pore Fluid Interactions; Chapman & Hall: New York, 1993. (4) Theng, B. K. G. The Chemistry of Clay-Organic Reactions; Wiley: New York, 1974. (5) Van Olphen, H. An Introduction to Clay Colloid Chemistry; Wiley: New York, 1963. (6) Sparks, D. L. Soil Physical Chemistry; CRC Press: Boca Raton, FL, 1988. (7) Sherwood, J. D. Langmuir 1994, 10, 2480-2486. (8) Delville, A.; Laszlo, P. New J. Chem. 1989, 13, 481-491. (9) Prezzi, M.; Monteiro, P. J. M.; Sposito, G. ACI Mater. J. 1997, 94, 10-17. (10) Kjellander, R.; Marcelja, S.; Pashley, R. M.; Quirk, J. P. J. Phys. Chem. 1988, 92, 6489-6492. (11) Norrish, K. Discuss. Faraday Soc. 1954, 18, 120-134.
unlimited swelling material, whereas calcium montmorillonite exhibits limited swelling due to ionic correlation forces.10-13 We have thus performed Monte Carlo simulations of the competitive condensation of counterions of different valency confined between two charged lamellae. Because of the deep electrostatic well in the vicinity of these uniformly charged surfaces, a major fraction of the neutralizing counterions are localized near the lamellae, leading to the so-called ionic condensation phenomenon. From the equilibrium ionic concentration profiles, we derived the ionic affinity of aqueous clay dispersions and the equation of state of clay suspensions. We determine also the minimum amount of monovalent counterions necessary to neutralize the cohesion of charged interface induced by divalent or trivalent counterions and to induce the crossover from cohesive to swelling behavior. In this study, the clay particles are modeled as smooth, infinite surfaces with uniform charge density. As shown previously,33 the spreading of ionic sites to form uniformly charged layers modifies the distribution of counterions only in the direct vicinity of the charged surfaces. The surface/ion and ion/ion interactions are described within the framework of the “primitive model”,14 characterized by long-range electrostatic couplings and contact repulsions. The main difficulty with numerical simulations performed within classical Euclidian geometry is the long range of Coulombian potential cut by the minimum image procedure. Ewald summation is the usual approach to introduce the long-range contribution of the electrostatic potential, but the accuracy of the results is a priori limited by the number of replica implied in the sum in the reciprocal space.15 Furthermore, self-consistent approaches using an external electric field16 are limited by the lateral extent of the ionic correlations occurring within (12) Delville, A.; Pellenq, R. J. M.; Caillol, J. M. J. Chem. Phys. 1997, 106, 7275-7285. (13) Pellenq, R. J. M.; Caillol, J. M.; Delville, A. J. Phys. Chem. B 1997, 101, 8584-8594. (14) Carley, D. D. J. Chem. Phys. 1967, 46, 3783-3788. (15) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford Science Publishers: Oxford, 1994. (16) Valleau, J. P.; Ivkov, R.; Torrie, G. M. J. Chem. Phys. 1991, 95, 520-532.
S0743-7463(98)00287-X CCC: $15.00 © 1998 American Chemical Society Published on Web 08/07/1998
5078 Langmuir, Vol. 14, No. 18, 1998
Figure 1. Scheme illustrating the use of hyperspherical geometry in order to model charged interfaces. The surface of the hypersphere accessible to the neutralizing counterions is limited by two parallel lamellae (labeled N and S) located at equal distances from the equatorial plane. The lamellae result from the intersections of the hypersphere by cones of angular aperture, θN and θS ) π - θN, respectively.
the layers of condensed counterions.13 To avoid these technical difficulties, we used a closed hyperspherical geometry,17 well suited to describe long-range electrostatic couplings. II. Methods The properties of charged interfaces are derived from equilibrium configurations of the neutralizing counterions (Figure 1). We study the competitive condensation of monovalent/ divalent and monovalent/trivalent counterions confined between two charged lamellae. Equilibrium configurations are generated by Monte Carlo simulations, using the classical Metropolis algorithm18 for sampling the (N,V,T) statistical ensemble. Ion-ion, clay-ion, and clay-clay interactions are described in the framework of the primitive model.14 The confined solvent molecules are considered as a continuum characterized by the dielectric constant of the bulk solvent (r ) 78.5 for water). In the present study, we have considered only the influence of ionic valency on the stability of charged interfaces. It should be interesting to investigate also the influence of the ionic radius, but all ions have here the diameter of solvated sodium (di ) 4.25 Å).19 Note, however, that the radii of solvated calcium (di ) 4.78 Å),19 magnesium (di ) 4.10 Å),19 and aluminum (di ) 3.8 Å)19 have the same order of magnitude. Such ionic size effects have been recently investigated.34 Clay particles are modeled as rigid and structureless infinite lamellae. Let us imagine the surface of a sphere drawn in a space of dimension 4. This 3D space is called hypersphere. Since the hypersphere is a closed space, it is possible to perform numerical simulations without any cutoff of the tail of the electrostatic potential. The hyperspherical geometry has been introduced to describe, in closed analytical form, the electrostatic couplings within homogeneous electrolyte modeled by charged hard spheres.17 We have recently applied the hyperspherical geometry to heterogeneous systems, to describe the interaction between charged lamellar interfaces12 (Figure 1), and predicted the stability of a large variety of homoionic materials.13 While the thermodynamic limit of the electrostatic energy is reached with a limited extension of the interface (less than 200 ionized sites), the derivation of the equation of state required some empirical correction in order to remove curvature effects.12,13 In the present study, we considered only two lamellar charge densities (σ ) 7 × 10-3 and 1.4 × 10-2 e/Å2) because they correspond to a large variety of minerals4,5 (montmorillonite, hectorite, vermiculite, mica, boehmite) and synthetic materials (laponite, metallic oxides [silicates, vanadium, or titanium oxides, aluminates]). These uniform charge densities would correspond (17) Caillol, J. M.; Levesque, D. J. Chem. Phys. 1991, 94, 597-607. (18) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087-1092. (19) Cooker, H. J. Phys. Chem. 1976, 80, 1084-2091.
Delville et al.
Figure 2. Calculated net force exerted across a fictitious plane separating the charged interface in two subunits (see text). to site-to-site separations varying between 8 and 12 Å if they were created by a regular cubic network. These separations are twice as large as the minimum interlamellar distance considered in the present study. The simulations were performed for 1000 elementary charges per lamella, with lateral extent in the interfaces varying between 260 and 380 Å. Starting from a random configuration of the confined counterions, 2.5 × 106 iterations are used to thermalize the ionic distributions, and 107 iterations are used to average the equilibrium properties of the interface. Typically, 24 h is required to perform one simulation on a workstation. The ratio of ionic contamination of the charged interface by monovalent counterions is defined by the average equivalent fraction of these monovalent ions:
x1° )
z1c1° z1c1° + z2c2°
(1)
where ci° is the average concentration of ion i (with valency zi) confined between the two charged lamellae. The same relationship is used to analyze the ionic composition of the diffuse layers by defining the local equivalent fraction of monovalent ions at two reference positions: either at direct contact with the charged interfaces (x1c) or at the plane midway between the two charged lamellae (x1f). These two positions are important since they quantify respectively the amount of condensed (x1c) and free (x1f) counterions. Obviously, the average ionic concentrations (ci°) are then replaced in eq 1 by the corresponding local ionic concentrations. Ionic activities (ai) and activity coefficients (γi) are approximated by the local ionic density at the plane midway between the lamellae:20
γi )
ai ci(D/2) ) ci ° ci °
(2)
The stability of the interface is derived from a local calculation of the longitudinal component of the pressure tensor:
P(z) ) (c1(z) + c2(z))RT +
1 S(z)
〈
∑∑ F R
Rβ‚ez〉
(3)
β
where ci(z) is the ionic density of ions i at position z, FRβ is the force exerted by particle β on particle R, and ez is the longitudinal director. The indexes (R,β) describe any pair of elements (ions or lamellae) located on different sides of the fictitious plane located at position z (Figure 2). This procedure was introduced by Berne and co-workers.21 S(z) is the interfacial area of the intersection between this reference plane and the hypersphere. The first contribution to the pressure (eq 3) is positive and results from the entropy of the interlamellar ions, and the second contribution (20) Marcus, R. A. J. Chem. Phys. 1955, 23, 1057-1068. (21) Rao, M.; Berne, B. J. Mol. Phys. 1979, 37, 455-461.
Stability and Ionic Affinity of Clay Particles
Langmuir, Vol. 14, No. 18, 1998 5079
is either attractive or repulsive12,17 and depends on the interparticular electrostatic and contact forces. Such derivation of the equation of state includes contributions from ionic correlations10,12,16 neglected by the classical Poisson-Boltzmann equation in the usual DLVO approach of colloidal stability. At equilibrium, the pressure is constant within the whole interface, and the calculated pressure is independent of the choice of the position z. By quenching the fictitious plane between one lamella and the ionic plasma,12 we obtain from eq 3:
P)-
σ2 + (c1(0) + c2(0))RT 20r
(4)
where ci(0) is the ionic density at contact with the lamella and σ is the charge density of the lamellae. The first contribution to the pressure is attractive and results from the electrostatic interactions between one isolated lamella and the other electric charges of the interface. The second contribution stems from the repulsive contact forces exerted by the condensed counterions on the lamella. No entropic term contributes explicitly to the pressure derivation. Because of the correlations between electrostatic and contact couplings, charged interfaces showed a great variety of attractive/repulsive behaviors,13 depending inter alia on the ionic diameter and valency, the interlamellar separation, and the surface charge density of the lamellae.
III. Results and Discussion (a) Competitive Condensation. Because of the strength of the electrostatic coupling, counterions are mainly localized in the vicinity of the charged lamellae (Figure 3a). This phenomenon is well known for polyelectrolytes and is called ionic condensation.22,23 However, by contrast with cylindrical24 or spherical25 polyions, no limiting law predicts the fraction of condensed counterions in the vicinity of lamellar polyions. Ionic condensation varies as a function of the interlamellar separation (D*, expressed in reduced units, is the interlamellar separation divided by the ionic diameter): At small separations (D* ) 1.18 in Figure 3a), ions spread uniformly within the available interlamellar space independently of their valency. At larger separations (D* ) 4.24 in Figure 3a), strong condensation is detected by differences of 1 order of magnitude or even more between the local ionic density at contact with the charged surfaces and at half-interlamellar separation. In the context of the central approximation which leads to the Poisson-Boltzmann equation and the DLVO theory, ions located at half-interlamellar separation are considered as free, and their local density is used to determine the ionic activity.20 Figure 3a also shows the large influence of the ionic charge on the strength of the condensation phenomenon. To better quantify the effect of electrostatic couplings on ionic condensation, we determine the equivalent fraction of contact (x1c) and free monovalent (x1f) counterions. As shown in Figure 3b,c, polyions exhibit a high preference for their most charged counterions: the increase of the electrostatic coupling (i.e., by increasing the ionic or lamellar charges) enhances the ionic selectivity of the interface. (b) Mechanical Stability. The equation of state of clay materials is drawn in Figure 4a-d. Electrostatic interactions are responsible for the attractive correlation (22) Manning, G. S. J. Chem. Phys. 1969, 51, 924-933. (23) Oosawa, F. Polyelectrolytes; M. Dekker Inc.: New York, 1971; Chapter 2. (24) Westra, S. W. T.; Leyte, J. C. Ber. Bunsen. Ges. Phys. Chem. 1979, 83, 672-677. (25) Delville, A.; Herwats, L.; Laszlo, P. New J. Chem. 1984, 8, 557562.
Figure 3. Analysis of the influence of ionic condensation, by a plot of (a, top) the ionic concentration profiles of monovalent (b, O) and trivalent (9, 0) counterions calculated for σ ) 0.014 e/Å2 and x1 ) 0.5; (b, middle) the local equivalent fraction of monovalent counterions at half-interlamellar separation (x1f, white symbols) and at contact with charged surface (x1c, black symbols) calculated with σ ) 0.007 e/Å2 at interlamellar separations D* ) 1.18 (circles), D* ) 2.25 (squares), and D* ) 4.24 (diamonds) for divalent competitive counterions; and (c, bottom) same as (b) but for σ ) 0.014 e/Å2.
forces between clay particles neutralized by divalent or trivalent counterions. At fixed interlamellar distance, the pressure continuously increases when divalent or trivalent counterions are gradually exchanged by monovalent
5080 Langmuir, Vol. 14, No. 18, 1998
Delville et al.
Figure 4. Equation of state of charged interface neutralized by mixed counterions at different fractions of monovalent cations (x1 ) 0 [O, s], 0.1 [3, - - -], 0.25 [4, - - -], 0.5 [0, - ‚ -], 0.75 [9, - - - - -], and 1 [b, ‚‚‚]). The lines are drawn to better visualize the results. Four interfacial systems were studied: (a, top left) surface charge density σ ) 0.007 e/Å2 with monovalent/divalent competitive cations; (b, top right) surface charge density σ ) 0.007 e/Å2 with monovalent/trivalent competitive cations; (c, bottom left) surface charge density σ ) 0.014 e/Å2 with monovalent/divalent competitive cations; (d, bottom right) surface charge density σ ) 0.014 e/Å2 with monovalent/trivalent competitive cations. Table 1. Limiting Fraction of Monovalent Counterions Necessary To Cancel Out the Cohesive Behavior of Charged Interfaces Neutralized by Polyvalent Cations surface charge density (σ), e/Å2
valency (z2)
limiting fraction (x1)
0.007 0.007 0.014 0.014
2 3 2 3
0.1 ( 0.05 0.5 ( 0.05 0.4 ( 0.05 0.5 ( 0.05
counterions. Figure 4 also defines the critical amount of monovalent ionic contaminant necessary to suppress the cohesive behavior of clay materials neutralized by divalent or trivalent cations (Table 1). As shown in Figure 5, the contribution from the monovalent counterions to the contact pressure (cf. eq 4) is monotonic and does not vary with the valency of the competiting counterions. The cohesive behavior of these charged materials results from the nonmonotonic variation of the contact pressure of divalent or trivalent counterions. The increase of the charge density of the lamellae emphasizes ionic condensation8 (cf. Figure 3) and the cohesive behavior of the charged materials (Figure 4 and Table 1). However, highly charged materials (σ ) 1.4 e/Å2 in Figure 4c,d) exhibit inflection points or local maxima of their equation of state if they are neutralized mainly by monovalent counterions. The location of these maxima (D* e 2) corresponds to the optimal maximum
Figure 5. Swelling pressure (9, 0) and its contact contributions from monovalent (2, 4) and polyvalent (1, 3) counterions. Black symbols correspond to divalent cations and white symbols to trivalent cations. The surface charge density σ ) 0.014 e/Å2, and the fraction of monovalent cation x1 ) 0.25.
of the interionic contact forces.13 To better illustrate this point, let us analyze the swelling pressure by locating the fictitious plane at half-interlamellar separation (z ) D/2 in eq 3). This procedure is more powerful since it allows
Stability and Ionic Affinity of Clay Particles
Figure 6. Swelling pressure (9) and its entropic (0), electrostatic (3), and contact (4) contributions calculated at halfinterlamellar separation for charge density σ ) 0.014 e/Å2 and fraction of monovalent cation x1 ) 1.
Langmuir, Vol. 14, No. 18, 1998 5081
Figure 8. Swelling variation of the chemical potential corresponding to the exchange of divalent (σ ) 0.007 e/Å2 [b]) or trivalent cations (σ ) 0.014 e/Å2 [9]) by equivalent amount of monovalent counterions (see text).
These results should quantify the efficiency of clay materials used as ionic exchangers.3 However, eq 2 neglects interionic correlations which contribute significantly to the stability of charged interfaces10,13,16 (see section IIIb). To better describe the chemical affinity of charged interfaces for monovalent counterions, we derive, from the equation of state, the variation of free energy characterizing the substitution of one polyvalent counterion by the equivalent amount of monovalent counterions. This ionic exchange mechanism maintains the overall electroneutrality of the interface. From eq 1 and classical thermodynamical relationships,26 we obtain in the (N,V,T) statistical ensemble: Figure 7. Swelling variation of the activity coefficients of monovalent (white symbols) and polyvalent (black symbols) counterions derived from eq 2. Circles and squares describe respectively interfaces with charge density σ ) 0.007 and 0.014 e/Å2.
us to distinguish the contributions12,16 due to electrostatic attraction and interionic contact repulsion. In the framework of the primitive model14 (charged hard spheres with finite radius), two sources of interionic correlations compete:12 electrostatic interactions, leading to attractive correlation forces, and contact repulsions, which are responsible for repulsive correlation forces (Figure 6). Coincidence between the positions of the local maximum of the swelling pressure and the maximum of the interionic contact repulsion is clearly shown in Figure 6. As a consequence, the increase of the surface charge density enhances two antagonistic contributions: electrostatic attraction and interionic contact repulsion. Highly charged materials have thus a more complex equation of state (cf. eqs 3 and 4), with coexistence of two phases.13 (c) Ionic Exchange and Affinity. The affinity of charged interfaces for monovalent ionic contaminants is another parameter characterizing the chemical stability of charged materials. Equation 2 is used to quantify the influence of competitive counterion condensation on the exchange capacity of the charged interfaces (Figure 7). Because of their stronger attraction within the electrostatic well in the vicinity of charged surfaces, polyvalent counterions have lower activity coefficients than monovalent counterions, except at small interlamellar separation, since ionic condensation then disappears (cf. Figure 3a-c). One also notes the influence of the surface charge density: doubling the interfacial charge increases the ion/ surface electrostatic coupling and the ionic condensation, thus reducing ionic activities.
µex(Vf) - µex(Vi) ) since
∂µex ∂P ∂2A )) ∂V ∂N1 ∂N1 ∂V
∂P dV ∫VV ∂N 1 f
i
(5a)
(5b)
In all cases, the pressure is a linear function of the equivalent fraction of the number of monovalent counterions (N1) at constant interlamellar separation. The slopes give access to ∂µex/∂V (cf. eq 5), which is integrated numerically. The result defines the chemical affinity of charged interfaces for exchanging polyvalent cations by monovalent counterions. This exchange free energy is not very sensitive to the valency of the polyvalent counterion (Figure 8). One further notes on Figure 8 a noticeable increase of the exchange free energy (∆µex ≈ 0.3-0.35) when the interlamellar separation is decreased from 1.6 to 1.2 reduced units. This result is not surprising, since the decrease of the interlamellar separation would not favor the exchange of polyvalent cations by monovalent because of the reduction of the available interlamellar space, which enhances excluded volume effects. (d) Swelling Energy. The variation of the energy of homoionic materials as a function of the interlamellar separation is shown in Figure 9. In all cases, swelling is an endothermic process, by contrast with the exothermic behavior measured for swelling clays (like sodium montmorillonite or laponite).27 Furthermore, the binding energy predicted for nonswelling clays (like sodium mica) is equal to -4.5 kJ/mol (see Figure 9), i.e., 1 order of (26) Hill, T. L. An Introduction to Statistical Thermodynamics; Dover Publishing Inc.: New York, 1960. (27) Fripiat, J.; Cases, J.; Franc¸ ois, M.; Letellier, M. J. Colloid Interface Sci. 1982, 89, 378-400.
5082 Langmuir, Vol. 14, No. 18, 1998
Delville et al.
are known to drastically alter the cohesive behavior of these materials after contamination by monovalent counterions. Obviously these predictions, based on the primitive model, are limited by the validity of this simple model. The primitive model neglects the atomic structure of the clay surfaces and the organization (preferential orientation,31 layering32) of the solvent molecules confined between the clay surfaces.29,30 As a consequence, this approach leads only to qualitative results for interlamellar separations lower than 4 solvent diameters (D* < 2.7). However, even in the framework of a molecular description of charged interfaces, a careful description of electrostatic interactions is still necessary, since the electrostatic couplings are then enhanced by the local reduction of the dielectric constant in the vicinity of the solid surfaces. IV. Conclusion Figure 9. Swelling variation of the energy of homoionic charged interfaces neutralized by monovalent (b, O), divalent (9, 0), and trivalent (4) counterions. Black and white symbols correspond respectively to surface charge density σ ) 0.007 and 0.014 e/Å2.
magnitude smaller than experimental data.28 These discrepancies of the primitive model have the same origin: the neglect of the atomic structure of the confined solvent and clay layers.29 The hydration energy of the clay surface and its neutralizing counterions30 contributes significantly to the swelling energy of clay materials and is responsible for clay swelling under the so-called “crystalline swelling regime”, i.e., for interlamellar separation smaller than 12 Å (D* e 2.7).11 (e) Limitations and Overview. The analysis of the equation of state of charged interfaces neutralized by mixed counterions clearly illustrates the influence of interionic contact repulsions (cf. Figure 6). As a consequence, this study should be extended by analyzing the stability of materials neutralized by counterions of different valency and diameter. More complex equations of state are expected to occur under such conditions, with multiple phase equilibria.13 One should also consider interfaces with higher charge densities (like cement surface under basic condition),13 for which alkali reactions9 (28) Gutshall, P. L.; Bryant, P. J.; Cole, G. M. Am. Miner. 1970, 55, 1432-1434. (29) Delville, A. Langmuir 1991, 7, 547-555. (30) Delville, A. J. Phys. Chem. 1993, 97, 9703-9712.
We used (N,V,T) Monte Carlo simulations in order to study the competitive condensation of monovalent/divalent and monovalent/trivalent counterions confined between two charged lamellae. The numerical simulations are performed within hyperspherical geometry in order to describe long-range electrostatic coupling in closed analytical form, without using any approximation. From the analysis of the ionic concentration profiles, the equation of state and the chemical affinity of charged interfaces were calculated by taking into account contributions from interionic correlations. We also derived the influence of monovalent ionic contaminant on the mechanical stability of charged interface. Acknowledgment. One of us N.G. thanks Dr. F. Bergaya and the “Fonds de la Coope´ration FrancoTunisienne” for financial support. The Monte Carlo simulations were performed either locally on workstations or on supercomputers (Cray C94 and C98 at IDRIS, Orsay, France). Our workstations were purchased thanks to grants from CNRS and Re´gion Centre (France). LA9802872 (31) Delville, A.; Grandjean, J.; Laszlo, P. J. Phys. Chem. 1991, 95, 1383-1392. (32) Israelachvili, J. N. Intermolecular and Surface Forces; Academic Press: London, 1985. (33) Soumpasis, D. J. Chem. Phys. 1978 69, 3190-3196. (34) Greberg, H.; Kjellander, R. J. Chem. Phys. 1998 108, 29402953.
