A Dissipative-Particle-Dynamics Model for Simulating Dynamics of Charged Colloid
aa r X i v : . [ c ond - m a t . s o f t ] N ov A Dissipative-Particle-Dynamics Model for Simulating Dynamicsof Charged Colloid
Jiajia Zhou ∗ and Friederike Schmid † Institut f¨ur Physik, Johannes Gutenberg-Universit¨at MainzStaudingerweg 7, D55099 Mainz, Germany
Abstract
A mesoscopic colloid model is developed in which a spherical colloid is represented by manyinteracting sites on its surface. The hydrodynamic interactions with thermal fluctuations aretaken accounts in full using Dissipative Particle Dynamics, and the electrostatic interactions aresimulated using Particle-Particle-Particle Mesh method. This new model is applied to investigatethe electrophoretic mobility of a charged colloid under an external electric field, and the influenceof salt concentration and colloid charge are systematically studied. The simulation results showgood agreement with predictions from the electrokinetic theory. ∗ [email protected]; † [email protected] . INTRODUCTION Colloidal dispersions have numerous applications in different fields such as chemistry, bi-ology, medicine, and engineering [1, 2]. In an aqueous solution, colloidal particles are oftencharged, either by ionization or dissociation of a surface group, or preferential adsorption ofions from the solution. A good understanding of the dynamics of charged colloids is impor-tant from the fundamental physics point of view. Furthermore, such an understanding mayalso provide insights to improve the material properties of colloidal dispersions. Theoreticstudies of charged colloids are difficult because of the complexity of the system and variousdifferent interactions among the colloid, solvents, and microions.Molecular simulations can shed light on the dynamic phenomena of charged colloids ina well-defined model system. Such studies are numerically challenging because two differ-ent types of long-range interactions are involved: the electrostatic and the hydrodynamicinteractions. In recent years, a number of coarse-grained simulation methods have beendeveloped to address this class of problem. The general idea is to couple the explicit chargeswith a mesoscopic model for Navier-Stokes fluids. One of the examples is the couplingscheme developed by Ahlrichs and D¨unweg [3], which combines a Lattice-Boltzmann (LB)approach for the fluid and a continuum Molecular Dynamics model for the polymer chains.This method has been applied to study the polyelectrolyte electrophoresis and successfullyexplained the maximum mobility in the oligomer range for flexible chains [4–6]. Besides theLattice-Boltzmann method [7–12], there are a few choices of the fluid model in the literature,such as the Direct Numerical Simulation (DNS) [13–16], Multi-Particle Collision Dynam-ics (MPCD) [17, 18], and Dissipative Particle Dynamics (DPD) [19–22]. In this paper, wechoose the DPD approach. DPD is a coarse-grained simulation method which is Galileaninvariant and conserves momentum. Since it is a particle-based method, microions can beintroduced in a straightforward manner. A recent comparative study [23] indicated thatthe electrostatic interaction is the most expensive part in terms of the computational cost.Therefore, for intermediate or high salt concentrations, different methods for modelling thefluid becomes comparable.For colloidal particles, one requirement of the simulation model is the realization of no-slip boundary condition on the colloid surface. In this work, we present such a colloid modelbased on the Dissipative Particle Dynamics, with full considerations of the hydrodynamic2nd electrostatic interactions. We apply this model to study the dynamics of a chargedcolloidal particle under static electric fields. The remainder of this article is organized asfollows: in section II, we introduce the simulation model and describe relevant parametersfor the system. We present the simulation results of electrophoretic mobility in section III.Finally, we conclude in section IV with a brief summary.
II. SIMULATION MODEL
Our simulation system consists of three parts: the solvent, the colloidal particle and themicroions. All simulations were carried out using the open source package ESPResSo [24].
A. Fluids
The fluids are simulated using Dissipative Particle Dynamics (DPD), an establishedmethod for mesoscale fluid simulations. For two fluid beads i and j , we denote their relativedisplacement as ~r ij = ~r i − ~r j , and their relative velocity ~v ij = ~v i − ~v j . The distance betweentwo beads is denoted by r ij = | ~r ij | and the unit vector is ˆ ~r ij = ~r ij /r ij . The basic DPDequations involve the pair interaction between fluid beads. The force exerted by bead j onbead i is given by ~F DPD ij = ~F D ij + ~F R ij . (1)The dissipative part ~F D ij is proportional to the relative velocity between two fluid beads, ~F D ij = − γ DPD ω D ( r ij )( ~v ij · ˆ ~r ij )ˆ ~r ij , (2)with a friction coefficient γ DPD . The weight function ω D ( r ij ) is a monotonically decreasingfunction of r ij , and vanishes at a given cutoff r DPDc , ω D ( r ) = (cid:18) − rr c (cid:19) if r ≤ r DPDc , r > r DPDc . (3)The cutoff radius r DPDc characterizes the finite range of the interaction.The random force ~F R ij has the form ~F R ij = q k B T γ
DPD ω D ( r ij ) ξ ij ˆ ~r ij , (4)3here ξ ij = ξ ji are symmetric, but otherwise uncorrelated random functions with zero meanand variance h ξ ij ( t ) ξ ij ( t ′ ) i = δ ( t − t ′ ) (here δ ( t ) is Dirac’s delta function). The fluctuation-dissipation theorem relates the magnitude of the stochastic contribution to the dissipativepart, to ensure the correct equilibrium statistics. The pair forces between two beads satisfyNewton’s third law, ~F ij = − ~F ji , hence the momentum is conserved. The momentum-conservation feature of DPD leads to the correct long-time hydrodynamic behavior (i.e.Navier-Stokes equations).In the following, physical quantities will be reported in a model unit system of σ (length), m (mass), ε (energy), e (charge) and a derived time unit τ = σ p m/ε . We use the fluiddensity ρ = 3 . σ − . The friction coefficient is γ DPD = 5 . m/τ and the cutoff for DPD is r DPDc = 1 . σ . To measure the shear viscosity, we implement the method in Ref. [25] tosimulate the Poiseuille and Couette flows in a thin channel geometry. The viscosity of thefluid with our parameter setting is η s = 1 . ± . m/ ( στ ). B. Microions
Microions (either counterions to balance the colloid charge or the dissolved electrolytes)are introduced as the same DPD beads as the fluid, but carry charges and have exclusiveinteractions (to other charged beads but not to the fluid beads). We only consider themonovalent case where microions carry a single elementary charge ± e . The exclusive inter-action is necessary to prevent the collapse of charged system. A modified, pure repulsiveLennard-Jones interaction is used [26], V ( r ) = ε "(cid:18) σr − r (cid:19) − (cid:18) σr − r (cid:19) + 14 if r − r ≤ r LJc , r − r > r LJc , (5)where r is the distance between two charged beads. The cutoff radius is set at the po-tential minimum r LJc = √ σ . The microions have a size of 1 . σ ( r = 0). Chargedmicroions also interact by Coulomb interactions, and we compute the electrostatic inter-actions using Particle-Particle-Particle Mesh (P3M) method [27–29]. The Bjerrum length l B = e / (4 πǫ m k B T ) of the fluid is set to 1 . σ and the temperature is k B T = 1 . ε .One useful quantity is the diffusion constant of the microion D I . We measure the mean-square displacement of the microion, then use a linear regression at late times to obtain the4iffusion constant, lim t →∞ h ( ~r ( t ) − ~r (0)) i = 6 D I t. (6)The diffusion constant depends on the salt concentration ρ s . We perform simulations withdifferent salt concentrations, ρ s = 0 . . σ − , in a simulation box L = 30 σ . Thesimulation results are compared with the empirical Kohlrausch law [30], which states thatmicroion’s diffusion constant depends linearly on the square root of the salt concentration √ ρ s , D I = A − B √ ρ s , (7)where A and B are two fitting parameters. Fig. 1 shows the simulation results and a fit toKohlrausch law. D I [ σ / τ ] ρ s [ σ -3 ]simulationfit 0.71-1.12 ρ s1/2 FIG. 1: The diffusion constant D I of microions as a function of salt concentration ρ s .Symbols are the simulation results. Error bars are obtained by averaging five independentsimulation runs. The curve is a fit to Kohlrausch law with fitting parameters A = 0 .
71 and B = 1 . C. Colloid
The colloidal particle is represented by a large sphere which has modified Lennard-Jonestype conservative interaction to the fluid beads. The interaction has a similar form of Eq.55), but with a larger radius R = r + σ = 3 . σ .To implement the boundary condition at the surface, a set of DPD interaction sites isdistributed evenly on the surface R = 3 . σ , and the number of the sites is N s . The positionof the interacting sites is fixed with respect to the colloid center. These sites interact withthe solvent beads through the DPD dissipative and stochastic interactions, with the frictionconstant γ CS and the same cutoff r DPDc . The total force exerted on the colloid is given bythe sum over all DPD interactions due to the surface sites, plus the conservative excludedvolume interaction, ~F C = N s X i =1 ~F DPD i ( ~r i ) + ~F LJ . (8)Here ~r i denotes the position of i -th surface sites. Similarly, the torque exerted on the colloidcan be written as ~T C = N s X i =1 ~F DPD i ( ~r i ) × ( ~r i − ~r cm ) , (9)where ~r cm is the position vector of the colloid’s center-of-mass. Note that the excludedvolume interaction does not contribute to the torque because the associated force pointstowards the colloid center. The total force and torque are then used to update the positionand velocity of the colloid in a time step using the Velocity-Verlet algorithm. The mass ofthe colloidal particle is M = 100 m and the moment of inertia is I = 360 mσ , correspondingto a uniformly distributed mass. Fig. 2 shows a representative snapshot of a single chargedcolloidal particle with counterions.As a benchmark for our colloid model, we have performed simulations of an unchargedcolloid in a simulation box L = 60 σ and measured the autocorrelation functions. Twofunctions are obtained from the simulations: the translational and rotational velocity auto-correlation functions C v ( t ) = h ~v (0) · ~v ( t ) ih ~v i , (10) C ω ( t ) = h ~ω (0) · ~ω ( t ) ih ~ω i , (11)where ~v ( t ) and ~ω ( t ) are the translational velocity and rotational velocity of the colloid attime t , respectively. Fig. 3 shows the simulation results for γ CS = 10 . m/τ in log-log plots.For short time lags, both autocorrelation functions show exponential relaxation. The6IG. 2: Snapshot of a colloidal particle in a salt-free solution. The surface sites arerepresented by the dark beads, and the light beads are counterions. For clarity, solventbeads are not shown here.decay rate can be calculated using the Enskog dense-gas kinetic theory [31–34]lim t → C v ( t ) = exp( − ζ v ENS t ) , (12)lim t → C ω ( t ) = exp( − ζ ω ENS t ) , (13)where the Enskog friction coefficients are ζ v ENS = 83 (cid:18) πk B T mMm + M (cid:19) / ρR M , (14) ζ ω ENS = 83 (cid:18) πk B T mMm + M (cid:19) / ρR M , (15)where m and M are the mass for the fluid bead and the colloid, respectively, and ρ is thesolvent density. Eqs. (12) and (13) are plotted as solid curves in Fig. 3 and show reasonableagreement with the simulation data when t < . τ .For long time lags, hydrodynamic effects set in and lead to a slow relaxation for auto-correlation functions [35]. This so-called long-time tail is the manifestation of momentumconservation, as the momentum must be transported away from the colloid in a diffusive7 -3 -2 -1 -2 -1 C v ( t ) t [ τ ]e -3.61t -3/2 -4 -3 -2 -1 -2 -1 C ω ( t ) t [ τ ]e -4.51t -5/2 FIG. 3: Translational (top) and rotational (bottom) velocity autocorrelation functions.The measurement is performed for a single uncharged colloid with radius R = 3 . σ in asalt-free solution. The temperature is k B T = 1 . ε and surface-fluid DPD parameter γ CS = 10 . m/τ .manner. Mode-coupling theory predicts an algebraic behavior at long times [36]lim t →∞ h ~v (0) · ~v ( t ) i = k B T mρ [ π ( ν + D )] / t − / , (16)lim t →∞ h ω (0) · ω ( t ) i = πk B Tmρ [4 π ( ν + D )] / t − / , (17)8here ν = η s /ρ is the kinematic viscosity and D is the diffusion constant of the colloid,which is much smaller than ν and can be neglected. The results are plotted as dashedlines in Fig. 3. The data are consistent with the theoretical prediction for t > τ , butthe rotational autocorrelation function exhibits large fluctuations for large times. This ismainly due to the fact that the statistics for long-time values becomes very bad, and verylong simulations are required to obtain accurate values.The diffusion constant of the colloid can be calculated from the velocity autocorrelationfunction using the Green-Kubo relation D = 13 Z ∞ d t h ~v (0) · ~v ( t ) i . (18)Alternatively, the diffusion constant can be obtained from the mean-square displacement,similar to Eq. (6). Due to the periodic boundary condition implemented in simulations, thediffusion constant for a colloid in a finite simulation box depends on the box size. Fig. 4demonstrates the finite-size effect by plotting the mean-square displacement as a functionof time for two different simulation boxes L = 10 σ and L = 30 σ . < [ r( t )-r( ) ] > [ σ ] t [ τ ]L=10 σ L=30 σ FIG. 4: The mean-square displacement of a spherical colloid with radius R = 3 . σ for twodifferent sizes of simulation box, L = 10 σ and L = 30 σ .The diffusion constant increases with increasing box size. For small simulation box, thelong-wavelength hydrodynamic modes are suppressed due to the coupling between the colloid9nd its periodic images. An analytic expression for the diffusion constant in terms of anexpansion of powers of 1 /L was derived by Hasimoto [37] D = k B T πη s (cid:18) R − . L + 4 . R L + · · · (cid:19) . (19)In Fig. 5, simulation results of the diffusion constant are plotted in terms of 1 /L , thereciprocal of the box size, and the curve is the prediction from Eq. (19). The simulationresults and the hydrodynamic theory agree well.
0 0.02 0.04 0.06 0.08 0.1 D [ σ / τ ] σ -1 ] Eq. (19)
FIG. 5: The diffusion constant D for a spherical colloid of radius R = 3 . σ as a function ofthe reciprocal of the box size 1 /L . The curve is the prediction from Eq. (19). Differentsymbols correspond to simulation runs with different initialization of the random generator.Recently, studies of flow over superhydrophobic surfaces demonstrate that no-slip bound-ary condition is not always appropriate [38, 39]. A more general boundary condition is theNavier boundary condition, where finite slip over the surface is allowed. One advantage ofour colloid model is the ability to adjust the boundary condition from no-slip to full-slip bychanging the surface-fluid DPD friction γ CS . Fig. 6 illustrates the change of the diffusionconstant by varying γ CS in a simulation box L = 30 σ . This freedom provides opportunitiesto study the effect of hydrodynamic slip on the colloidal dynamics [40, 41].10 .0080.0100.0120.0140.016 0.1 1 10 D [ σ / τ ] γ CS [m/ τ ] FIG. 6: The diffusion constant D for a spherical colloid of radius R = 3 . σ as a function ofthe surface-fluid DPD friction coefficient γ CS . The simulation box has a size of L = 30 σ . III. ELECTROPHORETIC MOBILITY
In this section, we apply our DPD-based colloid model to investigate the electrophoreticmobility of a charged colloid, and compare simulation results with predictions from elec-trokinetic theories [1, 42, 43].Charged colloids in an aqueous solution are surrounded by counterions and dissolved saltions. Counterions accumulate around the colloid surface due to the Coulomb attractionbetween opposite charges and form an electric double layer (EDL). In equilibrium, thecounterion cloud has spherical symmetry for a spherical colloid. The thickness of the EDLis characterized by the Debye screening length l D = κ − = " πl B X i z i ρ i ( ∞ ) − , (20)where the summation runs over different ion species; z i and ρ i ( ∞ ) are the valence and thebulk concentration for i -th ion, respectively. When an external electric field is applied to thesuspension, the colloid (assumed to be positively charged) starts to move in the direction ofthe field, while the counterion cloud (negatively charged) is deformed and elongated in theopposite direction of the field. A steady state is reached when the electric driving force isequal to the hydrodynamic friction acting on the colloid. When the field strength is small,11he final velocity of the colloid ~u depends linearly on the applied electric filed ~E , and theproportionality defines the electrophoretic mobility µ : ~u = µ ~E. (21)The electrophoretic mobility is in general a second-order tensor, but is reduced to a scalarfor spherical colloids. Due to the complexity of the system and the coupling between thehydrodynamic and electrostatic interactions, analytic solutions for the mobility only existfor limiting cases.In the literature, the electrophoretic mobility is often expressed in terms of a ζ -potential,defined as the electrostatic potential at the shear plane. In the limit of κR ≪ µ H = 2 ǫ m ζ η s , (22)where ǫ m is the fluid permittivity and η s is the shear viscosity. In the opposite limit whenthe Debye screen length is much thinner in comparison to the colloid size ( κR ≫ µ S = ǫ m ζη s = 32 µ H . (23)For more general cases of intermediate values of κR , one has to rely on numerical methodsto solve the electrokinetic equations [42, 43]. In this work, we compute the electrophoreticmobility using the software MPEK [54].In simulations, we use the colloid charge as a controlling parameter. To convert thecolloid charge to the zeta-potential, one need to solve the Poisson-Boltzmann equation. Fora spherical particle, numerical tables for the solution to Poisson-Boltzmann equation weregiven by Loeb et al. [46]. An analytic expression for the relationship between the ζ -potentialand the surface charge density σ was derived by Ohshima et al. [47, 48] σ = 2 ǫ m κk B Te sinh (cid:18) eζ k B T (cid:19) (cid:20) κR ( eζ / k B T )+ 1( κR ) eζ / k B T )]sinh ( eζ / k B T ) (cid:21) / . (24)We use the above equation to compute the zeta-potential from known surface charge density.12ig. 7 plots the electrophoretic mobility as a function of the dimensionless zeta-potential eζ /k B T . The simulations are performed for a charged spherical colloid of radius R = 3 . σ ,and the solution has a salt concentration ρ s = 0 . σ − . We follow the standard to use thereduced mobility, which is a dimensionless number and is defined as µ red = 6 πη s l B e µ. (25)The simulation results agree well with the prediction from the electrokinetic theory at smallzeta potentials. At large zeta potential, the steric effects of the microions may play a role[41, 49, 50], resulting in an increase of the mobility. µ r ed e ζ /k B Tsimulationmpek
FIG. 7: The reduced mobility as a function of ζ -potential for a spherical colloid of size R = 3 . σ and salt density ρ s = 0 . σ − .Fig. 8 plots the reduced electrophoretic mobility as a function of κR . In simulations, thecolloid has a fixed radius R = 3 . σ , and the salt concentration is varied to change the valueof κR . The charge on the colloid is small, Q = 10 e , to ensure that the zeta-potential isalso small. The simulation results and the theoretic predictions agree well, except at verylarge salt concentration (large κR value). One possible explanation to the discrepancy is thechange of microion’s diffusion constant with respect to the salt concentration (see Fig. 1),which is not taken into account in the electrokinetic theory.13 µ r ed κ R simulationmpek FIG. 8: The reduced mobility as a function of κR for a spherical colloid of size R = 3 . σ and charge Q = 10 e . IV. SUMMARY
We have developed a mesoscopic colloid model based on the Dissipative Particle Dy-namics. We have taken accounts in full for the hydrodynamic interaction with thermalfluctuations, using Dissipative Particle Dynamics, and the electrostatic interactions, usingParticle-Particle-Particle Mesh method. We applied this new colloid model to investigatethe electrophoretic mobility of a charged colloid under a static external field. The simulationresults show good agreement with the predictions from electrokinetic theories. Futhermore,this model has been applied successfully to study the dynamic and dielectric response of acharged colloid to alternating electric fields [51–53].
Acknowledgments
We are grateful to Prof. Reghan Hill for providing the computer program MPEK. Wethank the HLRS Stuttgart for a generous grant of computer time on HERMIT. This workis funded by the Deutsche Forschungsgemeinschaft (DFG) through the SFB-TR6 program14Physics of Colloidal Dispersions in External Fields”. [1] W. B. Russel, D. A. Saville, and W. Schowalter,
Colloidal Dispersions (Cambridge UniversityPress, Cambridge, 1989).[2] J. Dhont,
An Introduction to Dynamics of Colloids (Elsevier, Amsterdam, 1996).[3] P. Ahlrichs and B. D¨unweg, J. Chem. Phys. , 8225 (1999).[4] K. Grass, U. B¨ohme, U. Scheler, H. Cottet, and C. Holm, Phys. Rev. Lett. , 096104 (2008).[5] K. Grass and C. Holm, Soft Matter , 2079 (2009).[6] K. Grass and C. Holm, Faraday Discuss. , 57 (2010).[7] V. Lobaskin and B. D¨unweg, New Journal of Physics , 54 (2004).[8] V. Lobaskin, B. D¨unweg, and C. Holm, J. Phys.: Condens. Matter , S4063 (2004).[9] V. Lobaskin, B. D¨unweg, M. Medebach, T. Palberg, and C. Holm, Phys. Rev. Lett. , 176105(2007).[10] A. Chatterji and J. Horbach, J. Chem. Phys. , 184903 (2005).[11] A. Chatterji and J. Horbach, J. Chem. Phys. , 064907 (2007).[12] G. Giupponi and I. Pagonabarraga, Phys. Rev. Lett. , 248304 (2011).[13] H. Tanaka and T. Araki, Phys. Rev. Lett. , 1338 (2000).[14] Y. Nakayama and R. Yamamoto, Phys. Rev. E , 036707 (2005).[15] K. Kim, Y. Nakayama, and R. Yamamoto, Phys. Rev. Lett. , 208302 (2006).[16] Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E , 361 (2008).[17] A. Malevanets and R. Kapral, J. Chem. Phys. , 8605 (1999).[18] G. Gompper, T. Ihle, D. M. Kroll, and R. G. Winkler, Adv. Polym. Sci. , 1 (2009).[19] P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett. , 155 (1992).[20] J. M. V. A. Koelman and P. J. Hoogerbrugge, Europhys. Lett. , 363 (1993).[21] P. Espa˜nol and P. B. Warren, Europhys. Lett. , 191 (1995).[22] R. D. Groot and P. B. Warren, J. Chem. Phys. , 4423 (1997).[23] J. Smiatek, M. Sega, C. Holm, U. D. Schiller, and F. Schmid, J. Chem. Phys. , 244702(2009).[24] H. Limbach, A. Arnold, B. Mann, and C. Holm, Comp. Phys. Comm. , 704 (2006).[25] J. Smiatek, M. Allen, and F. Schmid, Eur. Phys. J. E , 115 (2008).
26] J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. , 5237 (1971).[27] R. Hockney and J. Eastwood, Computer Simulation Using Particles (Adam Hilger, Bristol,1988).[28] M. Deserno and C. Holm, J. Chem. Phys. , 7678 (1998).[29] M. Deserno and C. Holm, J. Chem. Phys. , 7694 (1998).[30] M. R. Wright,
An Introduction to Aqueous Electrolyte Solutions (Wiley, Chichester, 2007).[31] G. Subramanian and H. Davis, Phys. Rev. A , 1430 (1975).[32] J. Hynes, Annu. Rev. Phys. Chem. , 301 (1977).[33] J. T. Padding, A. Wysocki, H. L¨owen, and A. A. Louis, J. Phys.: Condens. Matter , S3393(2005).[34] J. K. Whitmer and E. Luijten, J. Phys.: Condens. Matter , 104106 (2010).[35] B. J. Alder and T. E. Wainwright, Phys. Rev. A , 18 (1970).[36] J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, London, 2006),3rd ed.[37] H. Hasimoto, J. Fluid Mech. , 317 (1959).[38] L. Bocquet and J.-L. Barrat, Soft Matter , 685 (2007).[39] O. I. Vinogradova and A. V. Belyaev, J. Phys.: Condens. Matter , 184104 (2011).[40] J. W. Swan and A. S. Khair, J. Fluid Mech. , 115 (2008).[41] A. S. Khair and T. M. Squires, Phys. Fluids , 042001 (2009).[42] R. W. O’Brien and L. R. White, J. Chem. Soc., Faraday Trans. 2 , 1607 (1978).[43] R. J. Hill, D. A. Saville, and W. B. Russel, J. Colloid Interface Sci. , 56 (2003).[44] E. H¨uckel, Phys. Z. , 204 (1924).[45] M. v. Smoluchowski, Z. Phys. Chem. , 129 (1917).[46] A. L. Loeb, J. T. G. Overbeek, and P. H. Wiersema, The Electrical Double Layer around aSpherical Colloid Particle (MIT Press, Massachusetts, 1961).[47] H. Ohshima, T. Healy, and L. White, J. Colloid Interface Sci. , 17 (1982).[48] H. Ohshima, Theory of Colloid and Interfacial Electric Phenomena (Academic Press, Ams-terdam, 2006).[49] J. L´opez-Garc´ıa, M. Aranda-Rasc´on, and J. Horno, J. Colloid Interface Sci. , 196 (2007).[50] J. L´opez-Garc´ıa, M. Aranda-Rasc´on, and J. Horno, J. Colloid Interface Sci. , 146 (2008).[51] J. Zhou and F. Schmid, J. Phys.: Condens. Matter , 464112 (2012).
52] J. Zhou and F. Schmid, Eur. Phys. J. E , 33 (2013).[53] J. Zhou, R. Schmitz, B. D¨unweg, and F. Schmid, J. Chem. Phys. , 024901 (2013).[54] http://reghanhill.research.mcgill.ca/research/mpek.html, 024901 (2013).[54] http://reghanhill.research.mcgill.ca/research/mpek.html