pubs.acs.org/Langmuir © 2010 American Chemical Society
Obtaining Effective Pair Potentials in Colloidal Monolayers Using a Thermodynamically Consistent Inversion Scheme A. D. Law* and D. M. A. Buzza Surfactant & Colloid Group, Department of Physics, University of Hull, Hull, HU6 7RX, United Kingdom Received October 1, 2009. Revised Manuscript Received January 31, 2010 The structure and stability of colloidal monolayers depends crucially on the effective pair interaction potential u(r) between colloidal particles. In this study, we construct a novel method for extracting u(r) from the two-dimensional (2D) radial distribution function g(r) of dense colloidal monolayers. The method is based on the Ornstein-Zernike relation and the HMSA closure first proposed by Zerah and Hansen (Zerah, G.; Hansen, J.-P. Self-consistent integral equations for fluid pair distribution functions: Another attempt. J. Chem. Phys. 1986, 84(4), 2336-2343). The HMSA closure contains a single fitting parameter which is determined by requiring thermodynamic consistency between the virial and compressibility equations of state. The accuracy of the HMSA inversion scheme is compared to a 2D predictor corrector scheme based on hard-disk fluids (HDPC) previously proposed by us (Law, A. D.; Buzza, D. M. A. Determination of interaction potentials of colloidal monolayers from the inversion of pair correlation functions: A two-dimensional predictor-corrector method. J. Chem. Phys. 2009, 131, 094704) and the conventional “one-step” inversion methods of HNC and Percus-Yevick (PY). The accuracy of all these schemes is tested against Monte Carlo simulation data for g(r) from monolayers interacting via a range of commonly encountered potentials, including both purely repulsive potentials and potentials containing an attractive well. For all the potentials studied, we find that the accuracy of the HMSA and HDPC schemes is superior to HNC and PY, especially as we go to higher densities. The HDPC and HMSA schemes are particularly accurate for hard-core and soft-core fluids, respectively, at high density and are therefore complementary to each other. Finally, we find that, even in the presence of experimentally realistic levels of noise in the input g(r) data, both HMSA and HDPC schemes are able to faithfully extract the salient features of the underlying interaction potentials. Both these schemes therefore provide convenient and accurate methods for extracting u(r) from experimental g(r) data for general 2D monolayers.
1. Introduction The properties of colloidal monolayers at interfaces have received significant attention in recent years due to their importance both scientifically, for example, in two-dimensional (2D) phase transition problems, and industrially, for example, in emulsion and foam stabilization, optical devices, and molecular electronic equipment.1,2 It is well-known that the structure and properties of colloidal monolayers are determined by the interaction forces between colloidal particles; understanding these forces is therefore vital for the applications above. However, unlike the case of three-dimensional (3D) colloidal suspensions, the interaction forces between colloidal particles at interfaces is very complex2 and still a matter of debate. In order to determine the actual nature of colloidal forces at the interface, it is critical that one is able to measure these forces accurately. A powerful method for obtaining the effective potential between particles is through the inversion of input structural data such as the radial distribution function, g(r), or the structure
factor, S(q).3-6 It should be emphasized that the potentials obtained from such inversions are effective pair potentials;7,8 that is, in addition to the bare pair potential, they may include many body effects (e.g., for dense systems) or external fields (e.g., for monolayers adsorbed on a substrate). Nevertheless, such effective potentials can be of considerable help in analyzing the behavior and properties of adsorbed monolayers. These inversion schemes make use of integral equation theory, where all nontrivial many body correlations are embodied within a function in the closure relation known as the bridge function, B(r).9 Unfortunately the bridge function is in general analytically intractable, and to make progress one therefore has to make approximations for B(r). The simplest inversion schemes are the so-called “one-step” inversion schemes where one assumes a simple form for the bridge function or, equivalently, a simple closure relation in the integral equation theory. The two most widely used closures are the Percus-Yevick (PY) and hypernetted chain (HNC) closures, and these have been used by a number of groups to invert pair-correlation data for 2D colloidal systems.3-5,10-12 However, it is widely recognized that
*To whom correspondence should be addressed. E-mail:
[email protected]. ac.uk. (1) Aveyard, R.; Clint, J. H.; Liquid droplets and solid particles at surfactant solution interfaces. J. Chem. Soc., Faraday Trans. 1995, 91, 2681-2697. (2) Bresme, F.; Ottel, M.; Nanoparticles at fluid interfaces. J. Phys.: Condens. Matter 2007, 19, 413101. (3) Brunner, M.; Bechinger, C.; Strepp, W.; Lobaskin, V.; von Grunberg, H. H. Density-dependent pair interactions in 2d. Europhys. Lett. 2002, 58(6), 926-965. (4) Quesada-Perez, M.; Moncho-Jorda, A.; Martı´ nez-Lopez, F.; Hidalgo-A lvarez, R. Probing interaction forces in colloidal monolayers: Inversion of structural data. J. Chem. Phys. 2001, 115 (23), 10897-10902. (5) Rajagopalan, R. Probing interaction forces in colloidal fluids through static structural data: The inverse problem. Langmuir 1992, 8, 2898-2906. (6) Rajagopalan, R.; Rao, K. S. Interaction forces in charged colloids inversion of static structure factors. Phys. Rev. E 1997, 55, 4423-4432.
Langmuir 2010, 26(10), 7107–7116
(7) Henderson, R. A uniqueness theorem for fluid pair correlation fluid. Phys. Lett. 1974, A49, 197-198. (8) Bolhuis, P. G.; Louis, A. A.; Hansen, J. P. Many-body interactions and correlations in coarse-grained descriptions of polymer solutions. Phys. Rev. E 2001, 64, 021801. (9) Hansen, J.-P.; McDonald, I. Theory of Simple Liquids.; Academic Press: London, 1986. (10) Likos, C. Effective interactions in soft condensed matter physics. Phys. Rep. 2001, 348, 267-439. (11) Baumgartl, J.; Arauz-Lara, J.; Bechinger, C. Like-charge attraction in confinement: myth of truth? Soft Matter 2006, 2, 631-635. (12) Behrens, S. H.; Grier, D. G. Pair interaction of charged colloidal spheres near a charged wall. Phys. Rev. E 2001, 64, 050401.
Published on Web 04/20/2010
DOI: 10.1021/la9037233
7107
Article
Law and Buzza
such closure relations are only approximate and may sometimes lead to inaccurate results. In a previous paper, we extended the hard sphere predictorcorrector method of Rajagopalan and Rao6 to two-dimensions and found that the resultant hard disk predictor-corrector (HDPC) scheme led to results that were much more accurate than the one-step inversion schemes for a wide range of potentials.13 However, one drawback of the scheme is that it is numerically unstable at densities close to the close-packed density of the reference hard-disk fluid. This is an important limitation in softcore colloidal monolayers where such densities are often encountered. To overcome this problem, in this paper, we propose an alternative inversion scheme based on thermodynamically consistent closure relations. The most well-known of these is the Rogers-Young (RY) closure14 which interpolates between the HNC and PY closures. The crossover from PY to HNC is governed by a switching function containing a single parameter, which is determined by imposing thermodynamic consistency. This “mixed” closure relation leads to very accurate results for the radial distribution function and thermodynamic properties in the case of purely repulsive potentials but breaks down for potentials containing an attractive well as thermodynamic consistency is no longer possible in this case.15 To overcome this problem, Zerah and Hansen15 have proposed an alternative mixed closure relation known as the HMSA closure which this interpolates between the HNC closure and the “soft-core” mean spherical approximation (SMSA) closure; the resultant closure is found to be applicable to both repulsive and attractive potentials. Zerah and Hansen have also used the HMSA closure to invert structural data from a 3D system interacting via the Lennard-Jones potential and achieved reasonable success. The aim of this paper is twofold. First, we extend the work of Zerah and Hansen and construct a HMSA based inversion scheme for two-dimensional monolayers. Second, we want to compare the accuracy of our HMSA scheme to the HDPC, RY, HNC, and PY schemes for a wide range of potentials. The rest of the paper is organized as follows. In section 2, we discuss details of the HMSA scheme. In section 3, we compare the accuracy of the HMSA to the HDPC, RY, HNC, and PY schemes by testing all these methods against g(r) data obtained from Monte Carlo simulations for a wide range of potentials (both repulsive and attractive). We also consider the effect of density, noise, and truncation of the r or q range of the input g(r) or S(q) data on the relative accuracy of these schemes; these are important issues to consider if we want to apply our inversion schemes to experimental data. Finally, in section 4, we present our conclusions.
2. HMSA Inversion Scheme The structural properties of an isotropic fluid are well-defined using a number of density correlation functions. An important example is the radial distribution function g(r), which is the average number density of particles at a radial distance r from any given particle divided by the average number density. The radial distribution function is related to the direct correlation function c(r) through the so-called Ornstein-Zernike equation9 Z hðrÞ ¼ cðrÞ þ F
c jr - r0 j hðr0 Þ dr0
ð1Þ
(13) Law, A. D.; Buzza, D. M. A. Determination of interaction potentials of colloidal monolayers from the inversion of pair correlation functions: A twodimensional predictor-corrector method. J. Chem. Phys. 2009, 131(9), 094704. (14) Rogers, F. J.; Young, D. A. New, thermodynamically consistent, integral equation for simple fluids. Phys. Rev. A 1984, 30, 999-1007. (15) Zerah, G.; Hansen, J.-P. Self-consistent integral equations for fluid pair distribution functions: Another attempt. J. Chem. Phys. 1986, 84(4), 2336-2343.
7108 DOI: 10.1021/la9037233
where h(r) = g(r) - 1. For actual calculations, it is more convenient to express the Ornstein-Zernike equation in Fourier space hðqÞ ¼ cðqÞ þ FhðqÞ cðqÞ
ð2Þ
where h(q) and c(q) are the Fourier transforms of h(r) and c(r), respectively. Equations 1 and 2 are exact, but since they introduce the new correlation function c(r), another relationship or closure is required in order to determine u(r) given g(r). Zerah and Hansen have proposed a “universal” closure relation for g(r) called the HMSA closure that interpolates between the HNC and the “soft-core” mean spherical approximation (SMSA) closures.15 The HMSA closure is universal in the sense that it applies to both repulsive potentials and potentials containing an attractive component. In the HMSA closure, we first separate the total interaction potential u(r) into a repulsive core part u1(r) and an attractive perturbative part u2(r) using the Weeks-ChandlerAnderson separation16 uðrÞ ¼ u1 ðrÞ þ u2 ðrÞ
ð3Þ
where ( u1 ðrÞ ¼
uðrÞ - uðrm Þ; r e rm 0; r g rm (
u2 ðrÞ ¼
uðrm Þ; r e rm uðrÞ; r g rm
and rm is the position of the primary potential minima, that is, the potential minima with the smallest r. The HMSA closure is then given in terms of u1(r) and u2(r) by g
HMSA
exp½f ðrÞ ðγðrÞ - βu2 ðrÞÞ - 1 ðrÞ ¼ exp½ - βu1 ðrÞ 1 þ ð4Þ f ðrÞ
where γðrÞ ¼ hðrÞ - cðrÞ
ð5Þ
f ðrÞ ¼ 1 - exp½ - Rr
ð6Þ
and
Here, f(r) is a simple switching function with the property f(r) = 0 for r f 0 and f(r) = 1 for r f ¥. The HMSA closure contains a single fitting parameter R in f(r) that is varied until thermodynamic consistency is achieved (thermodynamic consistency will be discussed in more detail later in this section). Note that eq 4 reduces to the Rogers-Young closure (RY)14 for purely repulsive potentials (i.e., u2 = 0 or rm = ¥). Having defined the HMSA closure, the inversion scheme based on this closure now proceeds as follows. The input data for the scheme is the radial distribution function g(r) (or the structure factor S(q)) and the corresponding number density F. In this paper, the input g(r) data is generated through Monte Carlo simulations (see next section), though our ultimate aim is to apply our inversion techniques to g(r) generated from experiment. (16) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54(12), 5237-5247.
Langmuir 2010, 26(10), 7107–7116
Law and Buzza
Article
In addition, we also need the correlation functions c(r), h(r), and γ(r) which are readily obtained from the input g(r) using the Ornstein-Zernike relation (eq 2). We use R = 0.1 as the initial guess for all cases studied. Equation 4 with the initial guess for R then provides the initial estimation of u(r). Specifically, for the region r g rm, u2(r) is given by eq 4 as ln½1 þ f ðrÞ hðrÞ - f ðrÞ hðrÞ ~ βu2 ðrÞ ¼ - CðrÞ - cðrÞ þ ð7Þ f ðrÞ Minimizing u2(r) with respect to r enables us to determine rm, the position of the primary minima of the potential. Having obtained rm, for r e rm, u1(r) is given by eq 4 as g~ðrÞ ~ mÞ - Cðr βu1 ðrÞ ¼ ln gðrÞ
ð8Þ
where
g~ðrÞ ¼
h i ~ m ÞÞ þ f ðrÞ - 1 exp f ðrÞðγðrÞ þ Cðr f ðrÞ
ð9Þ
Note that, in eqs7-9, we use h(r), γ(r), and c(r) calculated from the input g(r) data. Note also that eqs 7-9 correct for the typographical errors in the corresponding equations in the original Zerah and Hansen paper.15 Finally, the total potential is given by βuðrÞ ¼ βu1 ðrÞ þ βu2 ðrÞ
ð10Þ
An improved value of R is obtained by imposing thermodynamic consistency. Specifically, we require the isothermal compressibility calculated from the compressibility equation9 -1 Z ¥ DβP 1 ¼ ¼ 1 þ 2πF hðrÞr dr DF Sðq ¼ 0Þ 0
ð11Þ
to be equal to that calculated from the 2D virial equation17 βP πF ¼ 1F 2
Z
¥
r2 drgðrÞ
0
dβuðrÞ dr
ð12Þ
In order to calculate ∂βP/∂F via eq 11, we use the input h(r) on the right-hand side of eq 11. In order to calculate ∂βP/∂F via eq 12, we differentiate eq 12 with respect to F DβP ¼ 1 - πF DF
Z 0
¥
r2 drgðrÞ
dβuðrÞ πF2 dr 2
Z
¥
dr 0
DgðrÞ 2 dβuðrÞ r DF dr ð13Þ
The function ∂g(r)/∂F in eq 13 is then calculated by differentiating g(r) from eq 4 with respect to F to give 8 DγðrÞ > DgðrÞ < wðrÞ exp f ðrÞðγðrÞ - βuðrm ÞÞ DF ; r e rm ¼ DγðrÞ > DF : exp f ðrÞðγðrÞ - βuðrÞÞ ; r g rm DF
ð14Þ
where w(r) = exp(-βu(r)þβu(rm)). We now need to obtain ∂γ(r)/ ∂F, which can be readily obtained from the Ornstein-Zernike relationship (eq 2) and the HMSA closure (eq 4). In fact we have (17) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: Oxford, 1987.
Langmuir 2010, 26(10), 7107–7116
found that it is numerically more stable to work in terms of c(r) and γ(r), rather than c(r) and h(r). Rewriting eqs 2 and 4 in terms of c(r) and γ(r) and differentiating, we obtain DγðqÞ cðqÞ2 2FcðqÞ - F2 cðqÞ2 DcðqÞ ¼ 2 þ 2 DF DF 1 - FcðqÞ 1 - FcðqÞ
ð15Þ
and 8 > DγðrÞ DcðrÞ < DF wðrÞ exp f ðrÞðγðrÞ - βuðrm ÞÞ - 1 ; r e rm ¼ DγðrÞ > DF : exp f ðrÞðγðrÞ - βuðrÞ - 1Þ; r g rm DF ð16Þ Note that, in eqs 13-16, we have implicitly assumed that thermodynamic consistency is local, that is, that R (and hence f(r)) and u(r) are constant with respect to any changes in the density. The assumption that R is constant with respect to F is in general an excellent approximation, since R is always found to be a slowly varying function of density.15 Also note that, in eqs 15 and 16, we use the initial guess for R to calculate f(r) and the initial guess for u(r) calculated from eqs 7-10. We can now calculate ∂γ(r)/∂F by iterating between eqs 15 and 16 using a simple Picard scheme until convergence is obtained. For the initial guess of ∂γ(q)/∂F, we use the first term on the right-hand side of eq 15 only. The number of iterations required is approximately 200. This part of the scheme could be sped up by employing the Gillan-Newton-Raphson based method18 but we have found that the simpler Picard method is fast enough using a fast PC. To complete our calculation of ∂βP/∂F via eq 13, for g(r) in the second term of the right-hand side of eq 13, we insert the HMSA expression given by eq 4 with the value of R undetermined. An improved value of R is obtained by varying this R until ∂βP/∂F calculated from eq 11 agrees with ∂βP/∂F calculated from eq 13. This new value of R is now used to improve our estimate for u(r), and the procedure is repeated until R and u(r) converge. Typically, convergence of R to within 0.5% was achieved after about 5 iterations. All Fourier transforms were performed using the efficient method given by Lado,19 and typically the number of Fourier modes used was 1000. In order to avoid singularities in u(r) at r = 0, it is in fact numerically more convenient to work in terms of w(r) = exp(-βu(r) þ βu(rm)) instead of directly with βu(r) for the region r e rm. The relevant equations (i.e., eqs 8 and 13) in terms of w(r) rather than u(r) are given in Appendix A.
3. Results 3.1. HMSA versus One-Step Inversion Schemes. In this subsection, we compare the accuracy of the HMSA scheme with the HNC and PY schemes for a range of 2D potentials. We will also compare our HMSA scheme to the inversion scheme based on the RY closure, which is known to be very accurate for purely repulsive potentials. We will also compare the HMSA scheme to our previous hard-disk predictor corrector method (HDPC) in the next two subsections where we consider the effect of density, noise, and truncation on the accuracy of the different inversion schemes. Specifically, we consider the following four classes of potentials, which have been previously considered in ref 13: (18) Gillan, M. A new method of solving the liquid structure integral equations. Mol. Phys. 1979, 38, 1781. (19) Lado, F. Numerical fourier transforms in one, two and three dimensions for liquid state calculations. J. Comput. Phys. 1971, 8, 417-433.
DOI: 10.1021/la9037233
7109
Article
Law and Buzza
Figure 1. Inversion results for exponential decay potential, that is, eq A-3 with F* = Fσ2 = 0.015. Inset shows the comparison between the RY and HMSA inversions for this potential.
(1) (2)
(3) (4)
Simple exponential decay potential with and without an oscillatory tail; Stillinger-Hurd potential20,21 for colloids at interfaces which interpolates between a screened Coulomb potential at short-range and dipole-dipole interaction at long-range; Lennard-Jones (LJ) potential; Derjaguin-Landau-Verwey-Overbeek (DLVO) potential.
The full details of the potentials listed above can be found in Appendix B and in ref 13. There are a number of reasons for choosing these potentials. In addition to the fact that they are representative of the broad spectrum of interaction potentials found in 2D colloidal systems, they also contain additional attributes that allow us to test the inversion scheme rigorously. For example, the exponential decay potential and the StillingerHurd potential are purely repulsive and have a soft core. On the other hand, the exponential potential with oscillatory tail and the LJ potential have an attractive well, with the latter having a hard core. The DLVO potential is a very complex potential which includes a short-range hard core, an intermediate range attractive well, a longer range repulsive barrier, and finally a soft tail and therefore serves as a very stringent test of the inversion schemes. Finally, these potentials have been tested against our earlier HDPC scheme in ref 13 so that using these potentials allows for a direct comparison of the current HMSA inversion scheme with our previous HDPC scheme. As discussed in the previous section, we can use either g(r) or S(q) as the input data for our inversion scheme. Most experimental structural studies on 2D colloidal monolayers report results for g(r), and in this paper we therefore focus on g(r) as the input data. Following ref 13, g(r) data for the above potentials were generated using Monte Carlo (MC) simulations in the canonical ensemble using periodic boundary conditions and the minimum image convention. Simulations were performed with particle numbers ranging from 1024 to 2025 to confirm that finite size effects are negligible. For all systems studied, we chose densities high enough to exhibit several peaks in the g(r) data (20) Stillinger, F. H. Interfacial solutions of the poisson-boltzmann equation. J. Chem. Phys. 1961, 35, 1584. (21) Hurd, A. J. The electrostatic interaction between interfacial colloidal particles. J. Phys. A: Math. Gen. 1985, 18, L1055-L1060.
7110 DOI: 10.1021/la9037233
Figure 2. Inversion results for exponential decay with oscillatory tail potential, that is, eq A-4 with F* = Fσ2 = 0.015.
but low enough to ensure that the correlation length remains smaller than the simulation box size. In particular, the densities in this subsection were chosen to be low enough so that the first maximum in g(r) was significantly lower than 3, which is the value above which 2D fluids generally undergo an ordering transition to either a hexatic or crystalline phase;22 the specific densities used for each system are given below. However, in the next subsection where we consider the effect of density on the accuracy of the different inversion schemes, we will also consider higher densities close to ordering transitions for the monolayer. For all our simulations, 60 000 MC steps per particle were used for the equilibration phase while 100 000 MC steps per particle were used for the analysis phase. In this and the next subsection, g(r) curves are obtained by averaging over 10 000 snapshots from the latter phase of the simulation (i.e., 1 snapshot for every 10 MC steps per particle to ensure the snapshots are independent) in order to minimize noise. However, in subsection 3.3, where we consider the effect of noise in the input g(r) data on the accuracy of the inverted potentials, we will use g(r) data obtained from 1 snapshot to mimic experimentally realistic levels of noise in the input g(r) data. In collecting g(r) data, radial bin sizes of order 50 nm were used (assuming colloid diameters to be of order 1 μm); this level of resolution, though challenging, should be accessible experimentally. We first discuss results for the exponential decay potentials. Note that while the exponential decay potential is purely repulsive, the exponential decay with oscillatory tail has a primary minima around r/σ = 10 (see Figures 1 and 2). In both cases, a density of Fσ2 = F* = 0.015 was used, the same as that used by Quesada-Perez et al.4 Note that our previous HDPC method13 was numerically unstable at this density because the area fraction of the reference HD fluid was too close to close packing density and therefore a lower density had to be used. However, no such limitation was found for the HMSA scheme, as it is numerically stable even at this higher density. This indicates that the HMSA scheme is more suitable for soft potentials compared to the HDPC method because densities approaching the close packing density of the reference HD fluid are often encountered in these (22) In fact, for the 2D case, there is no universal value for the first maximum in g(r) where the fluid undergoes an ordering transition,23 that is, there is no 2D analog for the celebrated Hansen-Verlet criterion for crystallization in 3D.30 However, from available simulation data (e.g., refs 23 and 29), it would appear that 2D fluids generally enter an ordered phase when the height of the first maximum is around 3 or higher.
Langmuir 2010, 26(10), 7107–7116
Law and Buzza
Figure 3. Inversion results for screened coloumb and dipole potential, that is, eq A-5 with βζ = 1/5.75 10-4, ɛ = 80, and κ-1 = 200 nm and an area fraction φ = Fπσ2/4 = 0.005. Inset shows the comparison between the RY and HMSA inversions for this potential.
Figure 4. Inversion results for Lennard-Jones potential, that is,
eq A-6 with βζ = 1/1.25 and a reduced density of F* = Fσ2 = 0.5.
systems. This point will be discussed in more detail in the next subsection. The inversion of the exponential decay potential by all three schemes is shown in Figure 1. The HMSA scheme converges to a value of R = 0.130. Clearly, the HMSA scheme leads to a very accurate inversion of the potential and is much more accurate compared to either the HNC or PY. The inset of Figure 1 compares the inversion using the RY and HMSA methods. We see that the HMSA yields essentially the same result as RY. This is very encouraging, since the RY closure is known to be very accurate for purely repulsive potentials.10,14 In fact, strictly speaking, HMSA predicts a very shallow minima at rm = 28.462 with a depth βu(rm) = -0.016. However, since the magnitude of this minima is negligible compared to kBT, the inverted potential is effectively purely repulsive. In Figure 2, we show that the inversion of the exponential decays with an attractive well where the HMSA scheme converges to a value of R = 0.161. The scheme yields very accurate results as well, reproducing all the features, especially the primary minima and oscillatory tail, but slightly overestimating the soft core. In terms of the one-step methods, HNC overestimates the soft core Langmuir 2010, 26(10), 7107–7116
Article
Figure 5. Inversion results for DLVO potential, that is, eqs A-7-A-10, conducted at an area fraction φ = 0.2. Inset shows the source g(r) data used for this inversion.
to a much larger degree compared to HMSA, though it reproduces the primary minima and oscillatory tail well. In contrast, the inversion by PY is very poor for this potential, breaking down altogether between 8 < r/σ < 13. Thus, for potentials containing an attractive well, the potential is no longer bracketed by HNC and PY so it is not possible to use RY to invert the potential as thermodynamic consistency cannot be achieved in this case. In contrast, the HMSA scheme is able to faithfully invert potentials with an attractive well. Next, we examine the Stillinger-Hurd potential, which is a more realistic potential for colloids at a polar/nonpolar liquid interface. An area fraction φ = Fπσ2/4 = 0.005 was used. The inversion results are shown in Figure 3, and the HMSA scheme converges to a value of R = 0.170. Once again, the HMSA scheme leads to a very accurate inversion of the potential and is much more accurate compared to either the HNC or PY. We also see that the HMSA and RY schemes yield essentially the same inversion for this purely repulsive potential as shown in the inset of Figure 3. The third potential we test is the LJ potential, and the inversion results are shown in Figure 4. A reduced density of F* = Fσ2 = 0.5 was used for the monolayer. For the HMSA scheme, R converged to a value of R = 0.898. We see that the HMSA scheme reproduces the core and long-range tail accurately, although it overestimates the attractive well depth slightly. Both of the one step routines are also reasonably accurate at this density, though the HMSA scheme is marginally more accurate, especially in modeling the attractive well. The final potential we shall consider is the DLVO potential. As shown in Figure 5, this potential is very complicated as it includes a short-range hard core repulsion, an intermediate range attractive well, a longer range repulsion, and finally a soft tail and therefore serves as a very stringent test for any inversion scheme. The inversion was conducted at an area fraction of 0.2, and the results are shown in Figure 5. For the HMSA scheme, the converged value of R was R = 0.014. Note that this interaction potential required a higher number of Fourier modes for the inversion (approximately 2000 modes) because the sharp peak in the input g(r) data required higher resolution (see inset of Figure 5). We see that the accuracy of the HMSA scheme is reasonably good. Specifically, it predicts the short ranged hard core Born repulsion and the long-range soft tail well, though it it slightly underestimates attractive well depth and overestimates DOI: 10.1021/la9037233
7111
Article
Law and Buzza
Figure 6. Inversion results for exponential decay potential, that is, eq A-3 using all the considered inversion methods for (a) F* = Fσ2 = 0.008, (b) F* = 0.01, and (c) F* = 0.015.
the repulsive barrier, and its prediction of the position of the primary minima is slightly too large. For the one step routines, PY is more accurate compared to HNC, as it models the hard core repulsive, long-range tail, and primary attractive well more accurately. This is presumably due to the fact that PY is more accurate for short-range hard interactions compared to HNC and is therefore better adapted to the hard core repulsion present in the DLVO potential. However, PY predicts the attractive well depth less accurately compared to HMSA, though it predicts the repulsive barrier more accurately compared to HMSA. Taken as a whole, the accuracy of the HMSA scheme is therefore comparable to that of PY for this potential. In summary, for all potentials studied, we find that the HMSA scheme is at least as good as the most accurate one-step method for any given potential and, in most cases, is significantly more accurate. 3.2. Effect of Density on Accuracy of Inversion. In the previous subsection, we performed inversions at one specific density for each of the potentials studied. In order to study how the relative performance of the different inversion schemes depends on density, in this subsection, we compare the accuracy of the HMSA scheme with the HDPC, HNC, and PY schemes for a few representative potentials as we increase the density of the monolayer. We first consider the exponential decay potential without an oscillatory tail which has a soft core. In Figure 6, we present inversion results for this potential at the reduced densities of F* = Fσ2 = (a) 0.008, (b) 0.01, and (c) 0.015. We see that while HNC and PY yield reasonably accurate results at the lowest density, the accuracy of these schemes becomes progressively worse as we go to higher densities. On the other hand, both HMSA and HDPC yield very accurate inversion results for F* = 0.008 and 0.01 which are comparable to each other. However, as noted in the 7112 DOI: 10.1021/la9037233
previous subsection, the HDPC scheme breaks down at F*= 0.015 because, at this density, the area fraction of the reference hard-disk fluid (0.637 as determined from the HD equation of state13 and the compressibility equation (eq 11)) is too close to the close packed area fraction of η0 = 0.9069. In contrast, the HMSA scheme remains numerically stable at this density and still provides very accurate inversion results. Next, we consider the LJ potential which has a relatively hard core. In Figure 7, we present inversion results for the LJ potential at the reduced densities of F* = Fσ2 = (a) 0.4, (b) 0.6, and (c) 0.7. Just as for the case of the exponential decay potential, both HNC and PY yield accurate inversion results at the lowest density but become progressively worse as we go to higher densities. In contrast, the HDPC scheme yields accurate results for all the densities studied (though even this scheme becomes slightly less accurate at the highest density of F* = 0.7). Interestingly, while the HMSA scheme yields accurate results for F* = 0.4 and 0.6 (comparable to HDPC), it performs much worse than HDPC for F* = 0.7, though it is still marginally better than either HNC or PY. The poor performance of HMSA at F* = 0.7 is probably due to the fact that, at this density, the LJ fluid is close to an ordering transition. This is evidenced by the fact that the height of the first maximum in g(r) is about 2.7 which is close to the value of 3.3-3.4 where the LJ fluid undergoes an ordering transition to either a hexatic or crystalline phase.23 Since HMSA is fundamentally an integral equation theory for the isotropic fluid state, it is not surprising that the theory breaks down close to an ordering transition. In contrast, for the exponential decay potential, the height of the first maximum in g(r) is about 1.8 for the highest density studied. We therefore expect this system to be deep in the (23) Saija, F.; Prestipino, S.; Giaquinta, P. V. Entropy, correlations, and ordering in two dimensions. J. Chem. Phys. 2000, 113(7), 2806-2813.
Langmuir 2010, 26(10), 7107–7116
Law and Buzza
Article
Figure 7. Inversion results for Lennard-Jones potential, that is, eq A-6 for (a) F* = Fσ2 = 0.4, (b) F* = 0.6, and (c) F* = 0.7.
isotropic fluid state where the HMSA scheme is valid.22 Interestingly, the HDPC scheme yields accurate inversion results for all the densities studied for the LJ potential even though it is also based on integral equation theory.13 This suggests that the key approximation behind this scheme, that is, that the actual bridge function is equal to the hard disk bridge function, remains accurate even close to an ordering transition. Note that, unlike the exponential decay potential case, the HDPC scheme remains numerically stable at the highest density studied for the LJ potential. This is because the hard-core nature of the LJ potential drives long-range statistical ordering in the fluid. This means that, even close to an ordering transition, the area fraction of the reference HD fluid can still be quite far from the close packing density. For example, for F* = 0.7 where the height of the first maximum in g(r) is 2.7, the area fraction of the reference HD fluid is 0.605. In contrast, for the exponential decay potential, for F* = 0.01 where the height of the first maximum in g(r) is 1.52, the area fraction of the reference HD fluid is 0.552,24 while for F* = 0.015 where the height of the first maximum in g(r) is 1.71, the area fraction of the reference HD fluid is 0.637. Comparing our results for the exponential decay potentials and the LJ potential, we conclude that the accuracy of the HMSA and HDPC schemes are superior to HNC and PY, especially as we go to higher densities. For densities away from an ordering transition and the close packed density of the reference hard-disk fluid, both HMSA and HDPC schemes yield very accurate inversion results which are comparable to each other. For densities close to an ordering transition (as evidenced by the height of the first maximum in g(r) approaching the value of 3), the HDPC scheme is much more accurate compared to the HMSA scheme. The HDPC scheme is therefore better suited to studying hard-core (24) The entry for this quantity was erroneous in ref 13 due to a typographical error. We quote here the correct value.
Langmuir 2010, 26(10), 7107–7116
monolayers (e.g., monolayers interacting via the LJ potential) where there can be significant statistical ordering in the monolayer for relatively small area fractions of the corresponding hard-disk fluid. However, for densities close to the close packing density of the reference hard-disk fluid, the HDPC scheme becomes unstable while the HMSA scheme remains numerically stable. The HMSA scheme is therefore better suited to studying soft-core monolayers (e.g., monolayers interacting via the exponential potential) where significant ordering in the fluid only occurs for densities close to the close packing density of the reference fluid. The HMSA and HDPC schemes are therefore complementary to each other. 3.3. Effect of Noise and Truncation on Accuracy of Inversion. In the previous subsections, all inversions were performed on high quality g(r) data obtained by averaging over a large number of snapshots (10 000) in MC simulations. However, g(r) data obtained from real experiments will obviously be significantly more noisy. In order to quantitatively assess how errors in the input g(r) affect the final accuracy of the inversion, in this subsection we use the HMSA and HDPC schemes to invert more noisy g(r) data obtained from a single snapshot in our MC simulations. This procedure generates experimentally realistic levels of noise in the input g(r) data. In addition to noise, g(r) (S(q)) data obtained from experiments will in general be restricted in their r (q) range. We will therefore also consider how truncation in the r or q range affects the accuracy of the inversion. We first consider the exponential decay potentials, with and without an oscillatory tail, where the monolayers in both cases are at a density of F* = 0.015. Since the HDPC scheme is unstable at this density, we shall only consider inversions using the HMSA scheme for these potentials (we will consider the HDPC scheme later for the LJ potential). The advantage of considering these two potentials is that they have the same soft core but different intermediate and long-range behavior: the first has an intermediate DOI: 10.1021/la9037233
7113
Article
Figure 8. Radial distribution plots from simulations for exponential decay without an oscillatory tail (repulsive) and with oscillatory tail (attractive) potentials, both at F* = Fσ2 = 0.015. The solid and dashed lines represent g(r) data obtained by averaging over 10 000 snapshots for repulsive and attractive potentials, respectively, while the filled and unfilled circles represent g(r) data obtained from 1 snapshot for the repulsive and attractive potentials, respectively.
range attractive well, while the second has a purely repulsive tail. By comparing the two, we can therefore assess whether, in the presence of experimentally realistic levels of noise, our inversion scheme can still accurately distinguish between differences in the intermediate and long-range behavior of the underlying potential. It is particularly important that one is able to resolve differences in u(r) in this range, since there is considerable controversy in the literature over the existence of an intermediate range attractive component in the effective pair potential between quasi-2D charged colloids confined between two parallel plates.3,11,25-28 In Figure 8, we present g(r) data for these two potentials; the (solid and dashed) lines represent high quality g(r) data obtained from averaging 10 000 snapshots, while the data points represent the noisy g(r) data obtained from 1 snapshot. Clearly, the differences in g(r) between these two potentials is greater than the amplitude of the noise for the first peak and first trough, though the differences in g(r) become smaller than the noise amplitude for subsequent peaks and troughs. This shows that, even in the presence of noise, there is still sufficient information in the pair correlation function to distinguish between differences in the intermediate and long-range behavior of u(r). We also note from Figure 8 that the differences in the intermediate and longrange behavior of u(r) primarily manifest themselves as differences in g(r) in the intermediate r range, that is, in the first few maxima in g(r), while the differences at small and large r are small. This point is confirmed in Figure 9 where we plot (high quality) S(q) data for both the exponential decay potentials. In this case, the difference between S(q) primarily occurs in the intermediate q range, while the difference in the low q range (inset) is very small. The latter point is not surprising, since the low q regime of S(q) is controlled by the isothermal compressibility,9 which in turn is primarily determined by the core, rather than the long-range (25) Kepler, G. M.; Fraden, S. Attractive potential between confined colloids at low ionic strength. Phys. Rev. Lett. 1994, 73, 356-359. (26) Carbajal-Tinoco, M. D.; Castro-Roman, F.; Arauz-Lara, J. L. Static properties of confined colloidal suspensions. Phys. Rev. E 1996, 53, 3745-3749. (27) Crocker, J. C.; Grier, D. G. Microscopic measurement of the pair interaction potential of charge-stabilized colloid. Phys. Rev. Lett. 1994, 73, 352-355. (28) Crocker, J. C.; Grier, D. G. When like charges attract: The effects of geometrical confinement on long-range colloidal interactions. Phys. Rev. Lett. 1996, 77, 1897-1900.
7114 DOI: 10.1021/la9037233
Law and Buzza
Figure 9. High quality S(q) data for the exponential decay potentials without an oscillatory tail (solid line) and with an oscillatory tail (dashed) for F* = 0.015. The inset shows the results for the low q regime.
Figure 10. Inversion of exponential decay potential using the HMSA closure relation from g(r) data taken using only one snapshot.
region of u(r). From the above analysis, we conclude that while the core of u(r) influences the behavior of g(r) (S(q)) over the entire r (q) range, the influence of the intermediate and long-range region of u(r) on g(r) (S(q)) is primarily restricted to the intermediate r (q) range. Thus, while truncation of the r or q range of the input g(r) or S(q) data severely affects the accuracy of the inversion for the core region of u(r), we expect the effect to be much less severe for the intermediate and long-range region of u(r), provided that the truncated data includes the first few maxima in g(r) or S(q). In Figures 10 and 11, we present HMSA inversion results for the exponential decay potential without and with the oscillatory tail, respectively. In order to quantify the error in the inverted potential, for each potential, we inverted g(r) obtained from 3 independent snapshots. This resulted in 3 independent inverted potential outputs, from which we could calculate the mean and the standard deviation, which are the data points and error bars, respectively, shown in these figures. Clearly, despite the noisy input g(r) data, the HMSA inversion scheme is able to accurately distinguish between the purely repulsive tail in Figure 10 and the attractive well in Figure 11. This is particularly impressive given the very modest well depth of about 2kT used here for the Langmuir 2010, 26(10), 7107–7116
Law and Buzza
Article
Following the same procedure, in Figures 12 and 13, respectively, we present HMSA and HDPC inversion results for the LJ potential, where the monolayer is at a density of F* = Fσ2 = 0.5. Once again, despite the noisy input g(r) data, both the HMSA and HDPC schemes are able to accurately invert all the salient features of the LJ potential, including the repulsive core and the intermediate range attractive well. We note also that the accuracy of both schemes is comparable. The results in this subsection are very encouraging, demonstrating that both HMSA and HDPC schemes are robust with respect to noise and truncation of the input g(r) data and therefore provide a convenient and accurate method for extracting the effective pair interaction potential from experimental g(r) data for general 2D monolayers.
Figure 11. Inversion of exponential decay with oscillatory tail potential using the HMSA closure relation from g(r) data taken using only one snapshot.
Figure 12. Inversion of LJ potential using the HMSA closure relation from g(r) data taken using only one snapshot.
Figure 13. Inversion of LJ potential using the HDPC method from g(r) data taken using only one snapshot.
attractive potential. However, the amplitude of the noise in the input g(r) is such that it is not possible to resolve the sub-kT longrange oscillatory tail in u(r) using our HMSA scheme. Langmuir 2010, 26(10), 7107–7116
4. Conclusions In this paper, we have developed an inversion scheme for extracting effective pair potentials u(r) from 2D pair correlation functions g(r) which is based on the HMSA closure relation. In our implementation of the scheme, the unknown fitting parameter in the HMSA closure is determined by requiring thermodynamic consistency between the virial and compressibility equations of state. We have compared the accuracy of the HMSA scheme with our previous hard-disk fluid based predictor-corrector scheme (HDPC)13 and with inversion schemes based on hypernetted chain (HNC), Percus-Yevick (PY), and Rogers-Young (RY) closures for a range of 2D potentials, including exponential decay, Stillinger-Hurd, Lennard-Jones, and DLVO. We find that, for all these potentials, the HMSA and HDPC schemes are superior to HNC and PY, especially as we go to higher densities. In particular, the HMSA scheme is able to faithfully invert both purely repulsive potentials as well as potentials with an attractive well, and is therefore more universal compared to inversion schemes based on the RY closure. For densities close to an ordering transition, we find that HDPC scheme is more accurate than the HMSA scheme. The HDPC scheme is therefore better suited to studying hard-core monolayers where there can be significant statistical ordering in the monolayer for relatively small area fractions of the corresponding hard-disk fluid. On the other hand, the HMSA scheme remains numerically stable at densities close to the close packing density of the reference harddisk fluid where the HDPC scheme becomes unstable. The HMSA scheme is therefore better suited to studying soft-core monolayers where the monolayer can be far from an ordering transition even for densities close to the close packing density of the reference hard-disk fluid. In particular, the HMSA scheme is found to provide accurate inversion results provided the height of the first maximum in g(r) is significantly less than 3. The HMSA and HDPC schemes are therefore complementary to each other. Finally, we have considered the effect of noise and truncation of the r-range of the input g(r) data on the accuracy of the different inversion schemes. We find that truncation of the r-range of the input g(r) data affects the accuracy of the inversion much less for the intermediate to long-range region of u(r) compared to the core region of u(r). Thus, provided our primary focus is to resolve differences in u(r) in the intermediate to long-range (which for 2D colloidal monolayers it often is), we expect the HMSA and HDPC schemes to be reasonably robust with respect to truncation. We find that, even in the presence of experimentally realistic levels of noise in the input g(r) data, both the HMSA and HDPC schemes are able to faithfully extract the key salient features of the underlying interaction potentials. Both of these schemes therefore provide convenient and accurate methods for extracting the DOI: 10.1021/la9037233
7115
Article
Law and Buzza
effective pair interaction potential from experimental g(r) data for colloidal monolayers, and indeed 2D monolayers in general. Acknowledgment. The authors would like to thank the Surfactant and Colloid Group at the University of Hull for helpful discussions. A.D.L. also acknowledges the EPSRC for funding in the form of a DTA studentship.
Appendix A In order to avoid singularities in u(r) at r = 0 in our numerical calculations, it is more convenient to work with w(r) = exp(-βu(r)þβu(rm)) instead of directly with the potential βu(r) for the region r e rm. Rewriting the relevant equations in section 2.2 in terms of w(r) instead of βu1(r), eq 8 (for r e rm) becomes wðrÞ ¼
g~ðrÞ hðrÞ þ 1
ðA-1Þ
On the other hand, inserting eq 4 into eq 13 and using w(r) instead of βu1(r), eq 13 becomes Z rm DβP exp½f ðrÞðγðrÞ - βuðrm ÞÞ - 1 dwðrÞ ¼ 1 þ πF r2 dr 1 þ DF f ðrÞ dr 0 Z ¥ exp½f ðrÞðγðrÞ βuðrÞÞ 1 dβuðrÞ þ πF r2 dr 1 þ f ðrÞ dr rm 2 Z rm πF DγðrÞ dwðrÞ r2 dr exp½f ðrÞðγðrÞ - βuðrm ÞÞ þ DF dr 2 0 2Z ¥ πF DγðrÞ dβuðrÞ þ r2 dr exp½f ðrÞðγðrÞ - βuðrÞÞ ðA-2Þ DF dr 2 rm
The inversion now proceeds as described in section 2 but using the above equations in place of their corresponding equations in section 2.
Appendix B: Interaction Potential Details 1. Exponential Decay Potentials. Following QuesadaPerez et al.,4 the analytic form for the purely repulsive exponential decay is βuðrÞ ¼ 200expð- 0:65r=σÞ
ðA-3Þ
while the form for the exponential decay with oscillatory tail is βuðrÞ ¼ 200expð - 0:65r=σÞ þ 10expð- 0:15r=σÞ cosð0:3r=σÞ ðA-4Þ where σ controls the range of interaction.
7116 DOI: 10.1021/la9037233
2. Stillinger-Hurd Potential. For a colloid diameter of σ, we use the following potential
uðrÞ ¼
8
