Coarse-Graining Intermolecular Interactions in Dispersions of Highly

Letter pubs.acs.org/Langmuir

Coarse-Graining Intermolecular Interactions in Dispersions of Highly Charged Colloids Martin Turesson,*,† Bo Jönsson,‡ and Christophe Labbez*,† †

ICB, UMR 6303 CNRS, Université de Bourgogne, F-21078 Dijon Cedex, France Lund University, Theoretical Chemistry, POB 124 Chemical Center, 221 00 Lund, Sweden

‡

ABSTRACT: Effective pair potentials between charged colloids, obtained from Monte Carlo simulations of two single colloids in a closed cell at the primitive model level, are shown to reproduce accurately the structure of aqueous salt-free colloidal dispersions, as determined from full primitive model simulations by Linse et al. (Linse, P.; Lobaskin, V. Electrostatic Attraction and Phase Separation in Solutions of Like-Charged Colloidal Particles. Phys. Rev. Lett. 1999, 83, 4208). Excellent agreement is obtained even when ion−ion correlations are important and is in principle not limited to spherical particles, providing a potential route to coarse-grained colloidal interactions in more complex systems.

E

colloid charge densities. These effective parameters can be obtained by solving the Poisson−Boltzmann equation or by Monte Carlo (MC) simulations11−13 of a single spherical colloid in a Wigner−Seitz cell model or using the jellium approximation.14−17 However, at high electrostatic coupling, where ion−ion correlations are non-negligible, or when dealing with an anisotropic particle, this approach breaks down. In this letter, we show, on the basis of extensive MC simulations, that an accurate density-dependent wf(r) can actually be determined from sampling the mean force between two colloids in a cylindrical (closed) cell model (cyl-PM) at the PM level, supplemented with a Yukawa-like force at large separations. A few attempts were performed following a similar route but for very dilute suspensions in the weak coupling regime19 or with only qualitative agreement20 with the PM simulations. The coarse graining method introduced here is shown to give excellent results up to a reasonably high colloid volume fraction (∼7%) in both the low and high coupling regimes. Moreover, the method is not limited to isotropic particles. The ability to coarse grain the intercolloidal interactions is demonstrated (Figure 1) on the basis of a comparison between radial distribution functions g(r) obtained from the OCM simulations using the precalculated wf(r) and from the PM simulations of many-colloid systems with explicit monovalent and divalent counterions performed by Linse et al.18,21 In the cyl-PM and with the same system parameters as used in ref 18, the charged colloids were represented by impenetrable hard spheres with a diameter of σC = 40 Å carrying a central point charge of QC set to −60e, where e is the elementary charge. The simple ions were either monovalent, QI = e, or divalent, QI = 2e, with σI = 4 Å. The cyl-PM

lectrostatic interactions together with the shape, size, and concentration of colloidal particles are essential components in the design of self-assembling bulk materials with predefined properties. For example, electrostatic interactions have been shown to favor liquid crystal formation in dispersions of clays1 and gibbsite particles.2 The ability to tune electrostatic interactions makes it possible to form an ordered colloid crystal or a disordered gel phase from an isotropic phase and vice versa.3 However, although many beautiful experimental examples can be found in the literature,1−6 the modeling and understanding of such systems are still not mature. The main reasons are of course the system complexity and also the difficulty in performing direct computer simulations of a manyparticle system with charged colloids in the full primitive model (PM), which includes all ions explicitly. Already for small system sizes and moderate charge anisotropy, the sampling of the phase space becomes slow as a result of the large length scale difference between typical colloidal particles and simple ions.7 Because the simulation time scales as the square of the number of particles, an alternative route is to reduce the system to a one-component model (OCM), where only the colloids are present. The colloids, separated by a center-to-center separation r, interact via an effective pair potential wf(r) obtained through a thermal average over the simple ions. Such an approach has attracted wide interest and has been the focus of numerous studies. A well-known effective pair potential is the Yukawa potential8 acting between two charged colloids in a salt solution, derived from the Debye−Hückel theory. Despite its extensive use, the Yukawa potential, wY(r), is strictly valid only in dilute solutions, excess salt, and/or in the weak coupling regime. Inspired by the work of Oosawa and Manning9,10 on polyelectrolytes and the concept of “ion condensation”, the introduction of a densitydependent effective colloid charge and an effective Debye screening in wY(r) has been shown to give accurate wf(r) values for a wide range of colloid and salt concentrations as well as © 2012 American Chemical Society

Received: February 2, 2012 Revised: March 1, 2012 Published: March 9, 2012 4926

dx.doi.org/10.1021/la3005008 | Langmuir 2012, 28, 4926−4930

Langmuir

Letter

Table 1. Simulation Conditions in the cyl-PMa i

ϕi/%

L/Å

D/Å

Zf/e

κf−1/Å

1 2 3 4

0.84 1.7 3.3 6.7

500 400 300 250

142.73 112.84 92.13 71.36

[−18.5][−5.8] [−18.1][−5.4] [−18.5][−5.1] [−19.5][−5.4]

[61.4][99.8] [47.3][67.6] [34.1][42.3] [24.5][24.2]

a

The columns give (2) the colloid volume fraction, (3) the cylinder length, and (4) the cylinder diameter. The last two columns refer to the effective charges and screening lengths used in eq 4 for the Yukawa fits for monovalent and divalent ions, respectively.

Figure 1. For each volume fraction ϕ, the effective pair potential was determined from the sampled mean force (as a function of r) between two charged colloids fixed on the z axis in a cylindrical cell (middle image) in the presence of explicit counterions. For those ϕ and corresponding wf(r, ϕ) values, the radial distribution functions g*(r) in the OCM were subsequently simulated (right image) and compared with the exact g(r) obtained from the PM simulations18 (left image).

changes dramatically (Figure 2b) . The force profile now becomes nonmonotonic, with a deep attractive minimum at r ≈ σC + σI, followed by a decaying repulsive double-layer force at longer range that is around 10 times weaker in magnitude than for the system with monovalent ions. The attraction between like-charged objects, induced by the spatial correlation of ions on either surface, was first observed in the case of planar surfaces25,26 and later on in the case of macroions of different shapes.21,23,27−29 As in the case of monovalent ions, the consequence of increasing the volume fraction is a reduced magnitude of the repulsive barrier but also an increased depth of the minimum.30,31 A system size dependence was observed when sampling F(r, ϕ) at large r in the case of monovalent ions at the highest volume fraction. For divalent ions, however, this is not the case, showing that boundary effects are more pronounced when the ionic cloud around the colloids is more diffuse and thus more sensitive to the truncation at the cell boundaries. To access the long-ranged part of the intermolecular force, we fitted a function

simulations were performed at a temperature of T = 298 K with a uniform dielectric response of εr =78.4 with two charges i and j interacting, at a separation rij, by the Coulomb potential ⎧Q Q /4πε0εr rij , rij ≥ (σi + σj)/2 ⎪ i j u(rij) = ⎨ , rij < (σi + σj)/2 ⎪∞ ⎩

(1)

where ε0 is the permittivity of vacuum. The configurational space of the counterions was sampled by single particle displacements employing the Metropolis algorithm.22 The mean force at a given volume fraction F(r, ϕ) was evaluated over the midplane (z = 0) for a fixed colloid center-to-center separations r along the z axis. By doing so, the total mean force can be divided into three terms according to eq 2. F(r , ϕ) = F el(r , ϕ) + F id(r , ϕ) + F hs(r , ϕ)

F Y (r , ϕ ) = −

(2)

The electrostatic term Fel(r, ϕ), is calculated by summing all Coulomb forces between species residing on different sides of the midplane. The second term is the ideal contribution, which can conveniently be defined as βF id(r , ϕ) = [ρ I(z = 0) − ρ I(z = L /2)]A

exp[−κf r ] ⎛ 1 ⎞ ∂w Y (r ) ⎜ =A× + κf ⎟ ⎝r ⎠ ∂r r

(4)

derived from the Yukawa potential to the numerically calculated mean forces in the range of r = 75−100 Å (solid lines in the insets of Figure 2). In eq 4, A = lbZf2 exp[κfσC]/(1 + κfσC/2)2 is dependent on the effective colloid charge Zf and an effective screening length κf−1. Here, we stress the importance of the Yukawa fit to extrapolating the forces to large separations, which is required in the potential of mean force calculation and is given by r wf(r, ϕ) = −∫ ∞ F(r′, ϕ) dr′. As the volume fraction increases, κf−1 decreases because of the enhanced counterion concentration (screening) and the absolute value of Zf is significantly lower than the bare surface charge (Table 1). We now evaluate the accuracy of the current approach to calculating wf(r, ϕ) by comparing the radial distribution functions from MC simulations in the OCM, with those obtained in the PM simulations of the same systems.18 The OCM simulations were performed in a cubic simulation box using the minimum image convention with periodic boundary conditions applied in all three dimensions.32 Furthermore, the system was equilibrated through single colloid trial translations with the potential energy differences calculated to an arbitrary precision by using a cubic spline function between the points in the wf(r, ϕ) curve.33 The simulations were performed using 500 colloids, but no deviation in the calculated g(r) was found for a 10 times smaller system. The radial distribution functions with monovalent counterions in the OCM are shown in comparison with those of the PM in Figure 3a. Oscillating behavior is found

(3)

where ρI(z = 0) and ρI(z = L/2) are the ion densities at the midplane and at one cylinder wall and A is the cross-sectional cylinder area. The third term, Fhs(r, ϕ), stems from hard-core overlaps at the midplane and was calculated via a similar route suggested by Wu et al.23,24 The size dependence was checked for each system, but as a rule of thumb and besides the imposed colloid volume fraction (ϕ = 4/3 σC3D−2L−1) the cylinder length L was at least twice as large as the largest simulated colloid−colloid separation. To equilibrate the system, 5 × 104 MC trial displacements per ion were made and an additional 5 × 107 trial displacements to produce forces with absolute uncertainties of around 0.001 kT/ Å were made. The simulation parameters are given in Table 1. Figure 2 shows the calculated force curves in the presence of either monovalent or divalent counterions. In the investigated electrostatic coupling regime, we observe repulsive forces over the whole range of r in the case of monovalent counterions (Figure 2a). With increasing colloid volume fraction, the counterions are concomitantly concentrated, leading to an enhanced electrostatic screening that in turn lowers the doublelayer repulsion (inset). As the coupling strength is increased by replacing the monovalent ions with divalent ions, the picture 4927

dx.doi.org/10.1021/la3005008 | Langmuir 2012, 28, 4926−4930

Langmuir

Letter

Figure 2. Simulated mean forces (dashed lines with symbols) between two colloids at four different volume fractions as a function of r. (a) Monovalent counterions. (b) Divalent counterions. The insets show the corresponding Yukawa tail fits (solid lines). Note the different scales on the y axes.

Figure 3. Colloid−colloid radial distribution functions with monovalent counterions for various volume fractions (ϕ1, ϕ2, ϕ3, and ϕ4). (a) Solid lines are simulation results from the OCM using an effective potential, wf(r, ϕ), calculated in the cyl-PM. wf(r, ϕ) is obtained by a numerical integration of the force. For r ≤ 75 Å, the force is obtained by direct calculation (eq 2) and beyond by a Yukawa force fitted in the interval 75 ≤ r ≤ 100 Å (eq 4). Filled spheres are the corresponding PM data.18 (b) Solid lines as in part a,; the filled squares come from an OCM simulation with the Yukawa fit used over the whole range in r.

with the first peak appearing at a colloid−colloid separation of close to ρC−1/3, where ρC is the colloid number density. The agreement with the PM simulations18 is excellent for all but the highest colloid volume fraction, which displays a slight deviation. The observed discrepancy is due to the observed boundary effect in the cyl-PM, as previously mentioned. The positions of the local maxima and minima are, however, well reproduced for all volume fractions. In the case of monovalent counterions, the importance of the long-ranged part of the force curve is shown in Figure 3b, where radial distribution functions for the highest and lowest volume fractions have been obtained by extending the Yukawa force fits over the whole range of r prior to integration (solid lines in Figure 2a). The agreement with the g(r) obtained with wf(r, ϕ) from the cylPM is nearly perfect. This result clearly illustrates the reasons for the good accuracy previously obtained when using the effective Yukawa pair potential determined from a simple Wigner−Seitz cell model.12 Finally, the OCM radial distribution functions at higher electrostatic coupling with divalent counterions are displayed and compared to the exact g(r) in Figure 4a. Overall excellent agreement is found, demonstrating the high accuracy of the cyl-PM in predicting wf(r, ϕ), which implies that higher-order colloid correlations are of minor importance for the investigated volume fractions. As can be seen, the radial distribution functions are very different compared to the case of monovalent counterions, showing a growing peak with increasing ϕ at r ≈ σC + σI. This indicates the formation of colloid clusters, which were analyzed by

Figure 4. (a) Colloid−colloid radial distribution functions with divalent counterions for various colloid volume fractions (ϕ1, ϕ2, ϕ3, and ϕ4). Solid lines refer to the OCM, and filled circles refer to the PM simulation results.18 (b) Probability PN of finding clusters of colloids (simulation snapshots) with N members during an OCM simulation of the highest volume fraction with divalent counterions. Note the log scale on the y axis. 4928

dx.doi.org/10.1021/la3005008 | Langmuir 2012, 28, 4926−4930

Langmuir

Letter

defining a cluster as neighboring colloids with rij ≤ σC + σI (Figure 4b). For the highest volume fraction, clusters with up to six members were found, indicating microphase separation, in agreement with ref 21. In this case, approximating the potential of mean force over the full range of r with wY, as in Figure 3b, gives qualitatively wrong results. In summary, through the comparison between radial distribution functions of many charged colloids predicted from the OCM and the PM we have shown that the effective pair potential between colloids can be accurately determined from MC simulations, at the PM level, of two single colloids in a closed cylindrical cell. wf(r, ϕ) was obtained by integrating the sampled mean force versus the colloid separation, with the tail approximated by a Yukawa force at large separations. Overall, excellent agreement between the OCM and the PM is found when such a coarse graining of the colloidal interactions is employed. The agreement extends over a large range of volume fractions, up to ϕ ≈ 7%, but, more remarkably, it extends to the high coupling regime (i.e., when ion−ion correlations play a significant role). In view of the generality of the method, we expect this finding to be relevant to many complex systems going from dispersions and liquid crystals of charged anisotropic particles to the self-assembly of colloid particles at interfaces.

■

(7) Linse, P. Simulation of Charged Colloids in Solution. Adv. Polym. Sci. 2005, 185, 111. (8) Israelachvili, J. N. Intermolecular and Surface Forces; Academic Press: London, 1991. (9) Manning, G. S. Limiting Laws and Counterion Condensation in Polyelectrolyte Solutions I. Colligative Properties. J. Chem. Phys. 1969, 51, 924. (10) Oosawa, F. Polyelectrolytes; Marcel Dekker: New York, 1971. (11) Alexander, S.; Chaikin, P. M.; Grant, P.; Morales, G. J.; Pincus, P.; Hone, D. Charge Renormalization, Osmotic Pressure, and Bulk Modulus of Colloidal Crystals: Theory. J. Chem. Phys. 1984, 80, 5776. (12) Lobaskin, V.; Linse, P. Simulation of an Asymmetric Electrolyte With Charge Asymmetry 60:1 Using Hard-Sphere and Soft-Sphere Models. J. Chem. Phys. 1999, 111, 4300. (13) Denton, A. R. Poisson-Boltzmann Theory of Charged Colloids: Limits of the Cell Model for Salty Suspensions. J. Phys.: Condens. Matter 2010, 22, 364108. (14) Trizac, E.; Levin, Y. Renormalized Jellium Model for ChargeStabilized Colloidal Suspensions. Phys. Rev. E 2004, 69, 031403. (15) Castañ eda-Priego, R.; Rojas-Ochoa, L. F.; Lobaskin, V.; Mixteco-Sánchez, J. C. Macroion Correlation Effects in Electrostatic Screening and Thermodynamics of Highly Charged Colloids. Phys. Rev. E 2006, 74, 051408. (16) Dobnikar, J.; Castañeda-Priego, R.; von Grünberg, H. H.; Trizac, E. Testing the Relevance of Effective Interaction Potentials Between Highly-Charged Colloids in Suspension. New J. Phys. 2006, 8, 277. (17) Castañeda-Priego, R.; Lobaskin, V.; Mixteco-Sánchez, J. C.; Rojas-Ochoa, L. F.; Linse, P. On the Calculation of the Structure of Charge-Stabilized Colloidal Dispersions Using Density-Dependent Potentials. J. Phys.: Condens. Matter 2012, 24, 065102. (18) Linse, P.; Lobaskin, V. Electrostatic Attraction and Phase Separation in Solutions of Like-Charged Colloidal Particles. J. Chem. Phys. 2000, 112, 3917. (19) Allahyarov, E.; Löwen, H. Nonadditivity in the Effective Interactions of Binary Charged Colloidal Suspensions. J. Phys.: Condens. Matter 2009, 21, 424117. (20) Linse, P. Mean Force Between Like-Charged Macroions at High Electrostatic Coupling. J. Phys.: Condens. Matter 2002, 14, 13449. (21) Linse, P.; Lobaskin, V. Electrostatic Attraction and Phase Separation in Solutions of Like-Charged Colloidal Particles. Phys. Rev. Lett. 1999, 83, 4208. (22) Metropolis, N. A.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087. (23) Wu, J. Z.; Bratko, D.; Prausnitz, J. Interaction Between LikeCharged Colloidal Spheres in Electrolyte Solutions. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 15169. (24) Wu, J. Z.; Bratko, D.; Blanch, H. W.; Prausnitz, J. Monte Carlo Simulation for the Potential of Mean Force Between Ionic Colloids in Solutions of Asymmetric Salts. J. Chem. Phys. 1999, 111, 7084. (25) Guldbrand, L.; Jönsson, B.; Wennerström, H.; Linse, P. Electrical Double Layer Forces. A Monte Carlo Study. J. Chem. Phys. 1984, 80, 2221. (26) Kjellander, R.; Marčelja, S. Correlation and Image Charge Effects in Electric Double Layers. Chem. Phys. Lett. 1984, 112, 49. (27) Guldbrand, L.; Nilsson, L. G.; Nordensköld, L. A Monte Carlo Simulation Study of Electrostatic Forces Between Hexagonally Packed DNA Double Helices. J. Chem. Phys. 1986, 85, 6686. (28) Allahyarov, E.; D’Amico, I.; Löwen, H. Attraction Between LikeCharged Macroions by Coulomb Depletion. Phys. Rev. Lett. 1998, 81, 1334. (29) Woodward, C. E.; Jönsson, B.; Åkesson, T. The Ionic Correlation Contribution to the Free Energy of Spherical Double Layers. J. Chem. Phys. 1988, 89, 5145. (30) Grønbach-Jensen, N.; Beardmore, K. M.; Pincus, P. Interactions Between Charged Spheres in Divalent Counterion Solution. Physica A 1998, 261, 74.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]; christophe.labbez@ u-bourgogne.fr. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support from the Wenner-Gren Foundation and the support of the CRI from the University of Burgundy (https:// haydn2005.u-bourgogne.fr/CRI-CCUB) in accessing their computer facilities are gratefully acknowledged. We also thank Prof. Per Linse for providing us with the PM simulation data, together with Torbjörn Åkesson and Jan Forsman for valuable discussions.

■

REFERENCES

(1) Laurent, J. M.; Bihannic, I.; Maddi, S.; Funari, S. S.; Baravian, C.; Levitz, P.; Davidson, P. Liquid-Crystalline Aqueous Clay Suspensions. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 16101. (2) Mourad, M. C. D.; Byelov, D. V.; Petukhov, A. V.; de Winter, D. A. M.; Verkleij, A. J.; Lekkerkerker, H. N. W. Sol-Gel Transitions and Liquid Crystal Phase Transitions in Concentrated Aqueous Suspensions of Colloidal Gibbsite Platelets. J. Phys. Chem. B 2009, 113, 11604−11613. (3) Ruzicka, B.; Zaccarelli, E.; Zulian, L.; Angelini, R.; Sztucki, M.; Moussaid, A.; Narayanan, T.; Sciortini, F. Observation of Empty Liquids and Equilibrium Gels in a Colloidal Clay. Nat. Mater. 2011, 10, 56. (4) Kalsin, A. M.; Fialkowski, M.; Paszewski, M.; Smoukov, S. K.; Bishop, K. J. M.; Grzybowski, B. A. Electrostatic Self-Assembly of Binary Nanoparticle Crystals with a Diamond-Like Lattice. Science 2006, 312, 420. (5) Cölfen, H.; Antonietti, M. Mesocrystals: Inorganic Superstructures Made by Highly Parallel Crystallization and Controlled Alignment. Angew. Chem. 2005, 44, 5576. (6) Cölfen, H. Biomineralization: A Crystal-Clear View. Nat. Mater. 2010, 9, 960. 4929

dx.doi.org/10.1021/la3005008 | Langmuir 2012, 28, 4926−4930

Langmuir

Letter

(31) Lobaskin, V.; Lyubartsev, A.; Linse, P. Effective MacroionMacroion Potentials in Asymmetric Electrolytes. Phys. Rev. E 2001, 63, 020401. (32) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, 2002. (33) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in Fortran 77; Cambridge University Press: Cambridge, U.K., 1986.

4930

dx.doi.org/10.1021/la3005008 | Langmuir 2012, 28, 4926−4930