6718
Langmuir 1999, 15, 6718-6723
Simultaneous Monte Carlo Determination of Particle Size Distribution and Pair-Correlation Function of Spherical Colloids from a Diffraction Experiment Gergely To´th* Department of Theoretical Chemistry, Eo¨ tvo¨ s University Budapest, Budapest, P.O. Box 32, H-1518 Hungary Received February 12, 1999. In Final Form: May 24, 1999
A new simulation method was developed in which both positions and radii of the particles are treated as variables in a reverse Monte Carlo-like procedure. The new method can be used in the evaluation of diffraction data on one-component spherical colloids, where the polydisperse feature of the particles cannot be neglected. It reproduces the target functionsthe measured scattered intensitiesswell, without the a priori knowledge of the parameters of a supposed Schulz distribution of the particle size. The method is successful at a wide range of densities, but it fails at very dense systems and at broad distribution of particle size. Small-angle neutron and light-scattering data on spherical polystyrene colloids were evaluated by the new method and compared to data evaluated using the monodisperse assumption, where contact pairs were observed. A good reproduction of experimental data was achieved and no contact pairs were found in agreement with theoretical and other experimental results.
Introduction Wide-angle light scattering and small-angle neutron or X-ray scattering are efficient tools to get direct information on the structure of colloidal systems. There are two limitations to these methods: The accessible range in the inverse space is rather small, and the scattered intensity reflects both intra- and interparticle structure. The first limitation can be avoided by real-space models. In this case, the scattered intensity, I(q), of the model is calculated using Fourier transformation in the direction of the real to the inverse space. It avoids the error of Fourier transformation in the direction of inverse to real space containing truncation errors. The refinement of the model is based on the improvement of the reproduction of the measured function. There are several methods in the literature.1,2 One of them is the reverse Monte Carlo method3 that we have already applied for spherical colloids.4,5 The other problematic feature of the experimental data can be resolved if either the intra- or the interparticle contribution to the scattered intensity is neglected or approximated. Neglecting the interparticle interactions is rather common by measuring extra-dilute systems. In this case the interparticle structure is approximated with the ideal gas one, and the scattered intensity reflects solely the intraparticle structure. To evaluate the interparticle structure reflected in the structure factor, S(q), the scattered intensity should be corrected by the intraparticle part measured on extradilute systems or calculated on appropriate models. Unfortunately, the separation methods contain approximations; e. g., it is often assumed that the shape and the
size of the particles do not change with increasing concentration. Another feature of colloids biases the evaluation of diffraction data. Despite the more and more precise preparation of the particles, they never become identical as atoms or molecules. The particles have more or less well-defined size and shape distributions. The exact evaluation of data should be processed as on multicomponent systems, but the number of possible partials highly exceeds the number of independent measurements on the system despite the usage of contrast variation (isotopic substitution). It has been shown several times that inclusion of the polydispersity into the models improves the agreement of the calculated and the measured scattered intensities.6,7 Unfortunately, the size distribution should be built in the model a priori, and the trial and error refinement of the model becomes time-consuming. A mean structure factor can be evaluated to a mean paircorrelation function, g(r), with the exclusion of the polydispersity. In this case, the result becomes questionable, because it contains contact pairs4,5 that have not been observed by other experiments or model calculations. The scattered intensity depends strongly on the size distribution of the particles. In Figure 1, we show the scattered curves for an ideally dilute solution of polydisperse spheres (scatterers) for the case of a Schulz distribution with various half-widths (the functions are defined in the next section). Since the curves are characteristic with respect to the size distribution, we found it feasible to construct a reverse Monte Carlo-like procedure to optimize simultaneously the individual size and position of the particles in one simulation. Theoretical Background
* E-mail:
[email protected]. (1) Glatter, O. J. Appl. Crystallogr. 1977, 10, 415. (2) Glatter, O. In Neutron, X-ray and Light Scattering: Introduction to an Investigation for Colloid and Polymeric Systems; Lindner, P., Zemb, Th., Eds.; North-Holland: Amsterdam, 1991; p 33. (3) McGreevy, R. L.; Pusztai, L. Mol. Simul. 1988, 1, 369. (4) To´th, G.; Pusztai, L. J. Phys. Chem. 1992, 96, 7150. (5) To´th, G. Chem. Phys. Lett. 1997, 269, 413.
The normalized scattered intensity of a multicomponent system of particles with spherically symmetric scattering (6) Moonen, J.; de Kruif, C.; Vrij, A. Colloid Polym. Sci. 1988, 266, 1068. (7) D’Aguanno, B.; Klein, R. J. Chem. Soc., Faraday Trans. 1991, 87, 379.
10.1021/la9901534 CCC: $15.00 © 1999 American Chemical Society Published on Web 08/13/1999
Diffraction Experiment on Spherical Colloids
Langmuir, Vol. 15, No. 20, 1999 6719
f′i ) fi/∆Fsc and,
c ) cn(∆Fsc)2
Figure 1. Scattered intensities in arbitrary scale for ideally dilute systems (q* is in reduced units). 1: monodisperse spheres, 2-3: polydisperse spheres; 2: s ) 5%, 3: s ) 20%.
density and interparticle interactions is8 p
In(q) )
gij(r) ) g(r)
p
∑ ∑ fifjBi(q)Bj(q)(FiFj) i)1 j)1
1/2
Sij(q)
(1)
where fi is the scattering amplitude of the particles of component i, Bi(q) accounts for the intraparticle interference, p is the number of components, and Fi is the number density of component i. q ) 4πsin(θ/2)/λ, where λ is the wavelength and θ is the scattering angle. Using the approximation that the density of the scattering medium inside the particle is homogeneous and that the medium outside the particle equals that of the solvents, the scattering amplitude reads
fi ) 4π
∫0R ∆Fscr2dr
(2)
where ∆Fsc is the difference of the scattering medium density of the particle and the solvent and R is the radius of the i-th-type particle. Using this approximation, Bi(q) can be calculated analytically, because
Bi(q) ) f-1 i 4π
∫0R ∆Fscr2 sinqrqr dr ) 3(qR)-1[sin(qR) - qR sin(qR)] (3)
The Ashcroft-Langreth9 partial structure factors can be calculated from the partial pair-correlation functions, gij(r)-s:
Sij(q) ) δij + (FiFj)1/2
∫0∞ 4πr2(gij(r) - 1) sinqrqr dr
(4)
where δij is the Kronecker delta. The non-normalized intensity is a constant time the normalized one: p
I(q) ) cnIn(q) ) cn p
)c
which unifies the constants acting only on the amplitude of the scattered intensity. If we treat the system as an exact multicomponent one, the number of the possible partials is p(p-1)/2. In a polydisperse system, one may create fewer components by sorting the particles of similar sizes as the same type, using finite ∆R ranges. Unfortunately, in a numerical simulation for components at the edge of the size distribution, the number of particles is too small to give statistically satisfactory partial pair-correlation functions. On the other hand, the gij(r)-s of particles with relatively similar sizes do not differ significantly, at least in dilute solutions. For practical reasons we use the approximation:
p
fifjBi(q)Bj(q)(FiFj)1/2 Sij(q) ∑ ∑ i)1 j)1
p
f′i f′jBi(q)Bj(q)(FiFj)1/2 Sij(q) ) cI′(q) ∑ ∑ i)1 j)1
(5)
Here, (8) Fournet, G. In Small Angle Scattering of X-rays; Guinier, A., Fournet, G., Eds.; Wiley: New York, Chapman and Hall: London, 1955; p 65. (9) Ashcroft, N. W.; Langreth, D. C. Phys. Rev. 1967, 159, 500.
(6)
where g(r) is the pair-correlation function of treating the system as a one-component one. Of course, this approximation restricts the applicability of our method to systems where neither the polydisperse feature nor the packing fraction is too high. Method The reverse Monte Carlo method3 (RMC) is frequently used for the evaluation of diffraction data. The majority of its applications concern wide-angle scattering data on atomic and molecular systems,10 but there are also applications to colloids.4,5 The method uses the principles of the standard (Metropolis) Monte Carlo simulation, but instead of incorporating a priori known interparticle potentials, it moves the particles to reproduce a measured structural function, like S(q) or g(r). It does not apply the direct Fourier transformation of S(q) to g(r) to avoid truncation errors. The transformation is performed in the inverse direction, because the size of the simulation cell can be chosen large enough to reduce the truncation errors of the transformation. The simulation provides threedimensional configurations that can be analyzed for higher-order correlation functions, as long as the system can be described solely by pair interactions. The existence of these configurations also prohibits us from obtaining geometrically excludable partial pair-correlation functions, at least at the level of hard-sphere interactions. The reverse Monte Carlo method is one of the most effective tools available to derive pair-correlation functions from diffraction measurements. Furthermore, when the amount of the experimental data is not enough, constraints based on chemical intuition may be incorporated into the method (see the review of McGreevy).10 We modified the original algorithm to obtain an efficient tool for analyzing data on polydisperse colloids. We rewrote it to fit the configurations not to the structure factor, but to the scattered intensity, and we let the particle form factor change by varying the radii of the particles. Details of the algorithm are as follows: (1) Start with an initial configuration of N spheres in a periodic cube of sides L. The coordinates may be random or may correspond to some lattice. The radii are set to the experimentally suggested one. Calculate initial g(r), Bj(q), f′j(q), and I(q) for the configuration. Here j refers to (10) McGreevy, R. L. Nucl. Instrum. Methods Phys. Res., Sect. A 1995, 354, 1.
6720
Langmuir, Vol. 15, No. 20, 1999
To´ th
particles, since each particle is treated as a distinct component. The initial value of c is determined as, Nq
∑ i)1
c)
Nq
I′(qi) ∑ i)1
Iexp(qi)/
where Nq is the number of points in the measurement. Calculate Nq
∑ i)1
2
χ )
(
)
I(qi) - Iexp(qi) Iexp(qi)qi
2
(7)
where the denominator describes the approximate fluctuation of the scattered intensity originated from simulations with a finite number of particles.11 (2) Select one particle, move it by random displacement, and change its radius. Vary c randomly. Calculate g(r), Bj(q), f′j(q), I(q), and χnew2 for this new configuration. (3) Compare χnew2 and χ2. If χnew2 < χ2, the new configuration is accepted; otherwise, it is accepted with a probability proportional to Caccexp (χnew2-χ2). Cacc is adjusted to give a desired acceptance ratio. If the configuration is accepted, overwrite the value of χ2 by χnew2. (4) Repeat the procedure from (2) until χ2, c, and the size distribution of the particles have been converged. One can observe that the fiBi(q) weight of the i-th particle in the scattered intensity decreases as the radius of the particle decreases. If the radius of a particle decreases to a significantly smaller value than the radii of the other particles, the above algorithm does not effectively control the size and the position of this particle. Therefore, we have to build in an additional criterion. We calculated the distribution of the radii (hsim(r)) in NR discrete bins with the width of ∆R. We also calculated the mean (R) and the standard deviation (s) of the radii. Experimental evidence shows that the size distribution of polydisperse colloids can be well described by the Schulz distribution:
hideal(r) )
(t +R 1)
t+1
rt t+1 r exp R Γ(t + 1)
(
)
(8)
where t is the measure of the width and relates to the standard deviation of the distribution as
s ) (t + 1)-1/2R Γ denotes the gamma function. A χ2R value was calculated for each configuration that shows the deviation of the actual size distribution from an ideal Schulz distribution having the same mean and standard deviation: NR
χ2R )
∑ i)1
(
)
hsim (ri) - hideal (ri) hideal (ri)
2
(9)
The additional criterion was constructed similarly to the configuration described in step (3) and was applied simultaneously to the configurations generated in the algorithmic steps (1)-(4). The new configuration was not accepted if either the original or the additional criterion rejected it. We defined a separate caccR for the additional (11) To´th, G.; Baranyai, A. J. Chem. Phys. 1997, 107, 7402.
Figure 2. The success of the method in the test calculations with hard-sphere interactions. The s and the η values refer to the settings of the hard-sphere Monte Carlo simulations. 1: good reproduction, 2: weak reproduction, 3: false calculation. For explanation of the terms see the text.
criterion. We found that a 0.95 value for the partial acceptance ratio of the additional criterion prevented the undesired decreasing of the particle size and the loss of control of these particles. Test Calculations We had performed a test of the method on hard-sphere models before we applied it on real experimental data. Scattered intensities of various hard-sphere Monte Carlo simulations were calculated according to eqs 1-6. The simulation cells contained 1728 particles with fixed radii corresponding to Schulz distributions at given standard deviations. The number densities were calculated by approximating packing fractions (η) which were calculated using the mean of the distributions of the radii. Thereafter, we performed simulations with the new algorithm by fitting to the I(q)-s obtained by these hard-sphere Monte Carlo simulations. The number densities were chosen to be the same, the initial particle positions were chosen randomly, and the initial radii distribution corresponded to a monodisperse system with the same mean radius. The results of the test calculations are schematically shown in Figure 2. The reproduction of the scattered intensities (I(q)-s), pair-correlation functions (g(r)-s), distributions of the radii (h(r)-s), and such other functions as 〈Bi(q)〉 and S(q) were satisfactory in a reasonable range of η and s. We calculated the average deviation of these reproduced functions and the original functions. If the result was less than 2%, we denoted it as good in Figure 2. Between 2% and 20% we called weak, and over 20% we called a false reproduction. The method fails at either large η or s where the approximation gij(r) ) g(r) does not hold. The method also works well for broader distributions at medium densities, but it does not reproduce the data at low η and high s, despite that gij(r) ) g(r) holds here. It can be explained by the less-structured configurations of the hard-sphere model at this density that allow to reproduce the scattered intensities with more structured configurations, without reproducing the original size distribution. These simulations were unstable. Fortunately, in the case of the reproduction of real experiments detailed in the next section, the simulations could be performed without any unstable behavior, because the polystyrene particles can be modeled by long-range electrostatic interactions providing more structured configurations than the hard-sphere model at the same packing fraction. Therefore, we hope the method can be used at a larger range of η and s for systems where the hard-sphere interaction is not dominant.
Diffraction Experiment on Spherical Colloids
Figure 3. Particle-size distribution functions in arbitrary scale for the test calculations (reduced units are used: r* ) r/2R). 1-2: η ) 0.1, s ) 0.2. 1: new method, 2: input of the hardsphere Monte Carlo calculation, 3-4: η ) 0.05, s ) 0.05. 3: new method, 4: input of the hard-sphere Monte Carlo calculation.
Langmuir, Vol. 15, No. 20, 1999 6721
Figure 6. Particle-size distribution function in arbitrary units for polystyrene lattice obtained by the new method at volume fraction 0.13.
Figure 7. Scattered intensities in arbitrary units for polystyrene lattice at volume fraction 0.13. 1: fitted by the new method, 2: experimental data. Figure 4. Scattered intensities in arbitrary units for the test calculations (q* is in reduced units). 1-2: η ) 0.1, s ) 0.2.1: fitted by the new method, 2: calculated by hard-sphere MC method, 3-4: η ) 0.05, s ) 0.05. 3: fitted by the new method, 4: calculated by hard-sphere MC method.
at the beginning maximally with 5% and later with 1%. The maximal change in c was 1% in one step. A calculation took 2 days on a Pentium-II PC with 233 MHz. The Fortran code of the program is available on request to the author. Results on Experimental Data
Figure 5. Pair-correlation functions of test calculations (r* is in reduced units). 1-2: η ) 0.1, s ) 0.2. 1: new method, 2: hard-sphere MC method, 3-4: η ) 0.05, s ) 0.05. 3: new method, 4: hard sphere MC method.
The partial acceptance ratio for the distribution criterion was set to 0.95. The partial acceptance ratio of the main criterion was set to 0.1. We found that the dilute and less polydisperse systems converged slightly faster than the denser and more polydisperse ones. Five hundred accepted configuration/particle seemed to be enough to achieve convergence in all cases. In Figures 3, 4, and 5, we present the results of two simulations comparing to the original hard-sphere Monte Carlo data. The maximal random displacement was 0.17 times R, and the radii were changed
We selected two experimental data sets to reproduce by the new simulation method. Graphs of S(q) are more often published than I(q) for carefully investigated systems. We selected two measurements where S(q)-s had been published and the exact derivation of S(q)-s had been presented. We derived the scattered intensities from the structure factors. The S(q)-s were determined in both cases using the monodisperse approximation that provides the possibility of comparing the polydisperse evaluation of data to the monodisperse one of the traditional reverse Monte Carlo method. The first data set was obtained on a relatively dense (volume fraction of 0.13) polystyrene colloidal system by small angle neutron scattering measurement.12 The initial radii of the particles were set to the experimentally suggested 1.56 × 10-8 m for all of the 1728 particles. The functions of step (4) converged after a few hundred steps/ particle. The distribution of the radii shows the same mean as was proposed by the original article, but with s ) 8.3% (Figure 6). The reproduction of the experimental I(q) was satisfactory (Figure 7). We compared the g(r) to paircorrelation functions obtained by traditional reverse Monte Carlo simulations. In the latter cases, the system was treated as a monodisperse system (from the side of the form factors) and all the particles had the same radius. (12) Goodwin, J. W.; Ottewill, R. H. J. Chem. Soc. Faraday Trans. 1991, 87, 357.
6722
Langmuir, Vol. 15, No. 20, 1999
Figure 8. Pair-correlation functions for polystyrene lattice at volume fraction 0.13. 1: new method, 2-4: traditional reverse Monte Carlo simulation on a monodisperse system with a particle radius of 2: 1.56 × 10-8 m, 3: 0.94 × 10-8 m, 4: 0.312 × 10-8 m.
Figure 9. Particle-size distribution function in arbitrary units for polystyrene lattice obtained by the new method at volume fraction 0.0211.
As can be seen in Figure 8, the hard-core constraint had a strong effect on the results in the traditional reverse Monte Carlo. Using the experimentally suggested radius, we got a large peak for contact pairs which have an integral about 0.44 neighbor/particle. The peak was lowered and shifted to smaller distances using smaller hard-core radii. The hard-core radius effected the second and further peaks of the pair-correlation functions, too. This significant and unphysical dependence of data evaluation on the hardcore value shows the weakness of the monodisperse approximation of the traditional reverse Monte Carlo evaluation of data. In contrast, the pair-correlation function obtained by the new method does not imply contact pairs. Both the peak positions and the decay of the g(r) are similar, as can be observed by modeling with the Derjaguin-Landau-Verwey-Overbeek (DLVO) potential. The second experimental data set was measured by light scattering.13 It contained also polystyrene lattices at a lower volume fraction of 0.0211. The initial radii were set to the experimentally suggested 4.55 × 10-8 m, and hardsphere Monte Carlo random positions were applied at starting. The mean of the size distribution obtained by the new method was close to the originally suggested one (4.44 × 10-8 m) with s ) 6.2% (Figure 9). The reproduction of the scattered intensity was excellent (Figure 10). The large peak of the monodisperse data evaluation at small r disappeared if the new polydisperse method was applied (Figure 11). A small pre-peak still can be found in the g(r), (13) Tata, B. V. R.; Kesavamoorthy, R.; Sood, A. K. Mol. Phys. 1987, 61, 943.
To´ th
Figure 10. Scattered intensities in arbitrary units for polystyrene lattice at volume fraction 0.0211. 1: fitted by the new method, 2: experimental data.
Figure 11. Pair-correlation function for polystyrene lattice at volume fraction 0.0211. 1: new method, 2: traditional reverse Monte Carlo simulation on a monodisperse system with a particle radius of 4.55 × 10-8 m.
but we think it should be neglected. This pre-peak may be caused by experimental errors, multiple scattering, or inconsistent correction of data. Conclusion We presented a new numerical method in this paper that may serve as an efficient tool in the evaluation of diffraction data on spherical colloids. The most important feature is the simultaneous treatment of the intra- and the interparticle part of the experimental scattering intensities. In the presented case, we extended the reverse Monte Carlo procedure to the radii of the particles assuming a radius-dependent particle form factor. Of course, colloid particles are polydisperse not only in their size but also in their shape, e.g., spheres may distort to ellipsoids. This case is especially important for vesicles, but we did not consider this possibility in the present paper. Our results, obtained by the evaluation of experimental data, are not new. The polystyrene lattices can be satisfactorily modeled by DLVO potentials, and it is known that the polydisperse feature cannot be neglected. On the other hand, the previous models did not exactly reproduce the experimental scattered intensities. The results of previous reverse Monte Carlo calculations on monodisperse systems provided better reproduction of experimental data than the model calculations, but they showed the existence of contact pairs that were not proved by any other experimental or theoretical methods. Our results with the new numerical method are in good accordance with the theories and reproduce the experimental data excellently, without the questionable contact pairs.
Diffraction Experiment on Spherical Colloids
We should mention that it is not necessary to assume anything about the interaction potentials among the particles in the reverse Monte Carlo-like procedures. On the other hand, the reliability of the obtained configurations depends on the potentials. The main approximation of the present method is in eq 6. The method fails for systems with high density and polydispersity. Furthermore, simulations at small densities may become unstable, especially if the low packing fraction means enhanced
Langmuir, Vol. 15, No. 20, 1999 6723
configurational freedom. This is the case of hard-sphere systems, but it did not seem to be a danger if long-range (e.g., electrostatic) interactions are dominant. Acknowledgment. The author thanks the postdoctoral fellowship of the Hungarian National Science Foundation (OTKA). LA9901534
