J. Phys. Chem. 1987, 91, 219-222
219
Refined Monte Carlo Simulatlons of Static Percolationt J. Hoshen, R. Kopelman,* and J. S. Newhouse Department of Chemistry, The University of Michigan, Ann Arbor, Michigan 48109 (Received: June 6, 1986)
A modified Monte Carlo approach is used to derive and check the consistency of the critical site percolation concentrations and some critical exponents of the square and triangular lattices (with lo9 sites). While we confirm both universality and the accepted value for y and v (and c, for the triangular lattice), our c, value for the square lattice (c, = 0.5929 & 0.0002) is slightly higher than some recently quoted values. This reopens the, questions of methodology and lattice size. We also give effective y exponents for multiple layers, showing a monotonic crossover from two-dimensional to three-dimensional values. The preexponent is roughly unity over a wide concentration range.
Currently there is much interest in dynamic percolation problems and their relation to energy transport in molecular alloys, phonon (“fracton”) density of states, and heterogeneous chemical kinetics, giving fractal reaction orders for elementary The majority of simulations of fractal dynamics have been carried out on percolation cl~sters!~~However, only at exact criticality is the percolating cluster a fractal on all length scales. Thus, it turns out to be of much practical importance to know the exact value of the critical concentration (c,) of square (and also cubic) lattices. We note that even for the triangular lattice, we cannot assume that the theoretical value c, = 1/2 is valid for finite (though large) lattices used in the Monte Carlo simulations. Furthermore, a serious question arises whether different criteria result in the same c,, especially for finite lattices. The Monte Carlo method has been extensively used in percolation studies since this phenomenon has been described by Broadbent and Hammersley in 1957.6 This technique has been applied extensively to determine a variety of percolation parameters, such as exponents, cluster sizes, and cluster shape parameters. A comprehensive, up-to-date summary of the field is provided in the just appeared monograph by S t a ~ f f e r . The ~ main advantage of the Monte Carlo method is that it gives results when other methods are not available. Yet, there are some lingering questions. Namely, can the simulation results be accepted as quantitative measures for the percolation parameters? In the next section we shall try to provide answers to this question, drawing heavily on our own experience. One of our aims has been to investigate, and possibly improve, the c, value of the two-dimensional lattices, using a somewhat new approach, Le., utilizing the new theoretical values7 for the twodimensional static critical exponents (@= 5/36, y = 43/18, etc.). Conversely, having established good values for the finite lattice critical concentration, we test the validity of the scaling relation for finite lattices. In particular we derive the exponents y and v. Furthermore, we investigated the dependence of the critical exponent y on the number of two-dimensional layers (moving, quasi-continuously, from the two-dimensional square lattice to the three-dimensional cubic lattice). This is also of some practical interest,* as is the value of the preexponent.
Monte Carlo Method We shall briefly review the Monte Carlo method for the percolation site problem, where sites are occupied with a probability cg (or are unoccupied with a probability of 1 - cg). The Monte Carlo method uses random numbers to simulate a random lattice, where a site n is marked occupied if a random number X,, such that 0 < X , < 1, is less than cg. In our simulation, we define a random lattice configuration as the set of lattice sites generated by an identical random number vr,l sequence derived from a single initial seed. For a given configuration, new occupied sites are created, with the increase of cg (Occupied sites previously created at lower cg do not disappear or change their location.) Enumerating the number of clusters and their sizes turns out to be +This work was supported by NSF Grant DMR 8303919.
0022-3654/87/2091-0219$01.50/0
TABLE I: Y for Triangular Monolaver Lattices edge length Y correlation ~
CC
4 000 16 000 32 000
0.5001 0.5001 0.5002 0.5001 0.5000
2.25 2.331 2.457 2.404 2.351
0.999 99 0.99994 0.999 987 0.999 996 0.999999
configurations 5 1 1
1 1
the more complicated part of the simulation. There are a variety of methods for the enumeration.’ We have been using the CMLT7q9*10 method, which proved to be exceedingly efficient for very large simulated lattices, comprising billions of sites. The Monte Carlo run gives the cluster size distribution for the simulated lattice. In other words, we obtain the number of clusters i , in the lattice for every cluster size m . From the distribution of the cluster sizes i,, we can evaluate the average cluster size, raw
1 I,” = -Ei,m2 G where the summation is taken over all cluster sizes. Another useful quantitys-10 is Iayl which is the average cluster size, excluding the contribution of the largest cluster, “ax
2
I,,,’ = I,, - G where m,,, is the size of the largest cluster (maxi cluster). I,,,’ rises very rapidly below the critical percolation point c, and drops even faster above the critical point. We use the maximum of the curve Ia,,’,plotted as a function of cg, to determine c,.
Cluster Size Exponent y According to the scaling theory,’ the average cluster size I,, below the critical point is given by cg - c,
Iav
= k(
-’
7)
(3)
where y is a critical exponent and k is a constant. A log-log plot of (3) will give a straight line, where y is the line slope. Obtaining y seems to be rather a simple procedure, yet getting an accurate (1) Kopelman, R. Modern Problems in Solid State Physics; Agranovich, V. M., Hochstrasser, R. M., Eds.; North Holland: Amsterdam, 1983;vol. 4, p 139. (2) Alexander, S.;Orbach, R. J . Phys. Lett 1982, 43, L625. (3) Kopelman, R.; Parus, S.; Prasad, J. Phys. Reu. Lett. 1986, 56, 1742. (4) Argyrakis, P.; Kopelman, R. Phys. Reu. B 1984, 29, 51 1. (5) Keramiotis, A.; Argyrakis, P.; Kopelman, R. Phys. Rev. B 1985, 31, 4617. (6) Broadbent, R.; Hammersley, J. M. Proc. Cambridge Philos. SOC.1957, 53, 629. (7) Stauffer, D. Introduction to Percolation Theory; Taylor & Francis: New York, 1985. (8) Hoshen, J.; Kopelman, R.; Monberg, E. M. J. Stat. Phys. 1978, 19, 219. (9) Hoshen, J . ; Kopelman, R. Phys. Rev. B 1976, 14, 3438. (10) Hoshen, J.; Klymko, W.; Kopelman, R. J. Stat. Phys. 1979, 21, 583.
0 1987 American Chemical Society
The Journal of Physical Chemistry, Vol. 91, No. I, 1987
220
Hoshen et al.
B -
7 6 -
b
w
5 -
: m 4 -I
3
2 1
0
0
1
2
-Log
I
(CC
3
- Cg)
a
1
0
2
4
3
/ ccl
g 1
-Log
I
(CC
- Cg)
/ Ccl
-Log
1 (Cc
- Cg)
/ ccl
Figure 1. Derivation of y for a square lattice. A log-log plot for the average cluster size vs. I(c - c,)/c,l. Lattice sizes: (a, top left) 400 (b, top right) 4000 X 4000 sites, (c, bottom left) 16000 X 16000 sites, and (d, bottom right) 32000 X 32000 sites.
of (3) turns out to be a formidable task. The two major problems are that y is very sensitive to even small changes in c, and near the critical point the values of I,, exhibit a significant scatter. In Figure'l we have a log-log plot of Z, vs. (cg - c,)/c, for lattices extending from 400 X 400 to 32000 X 32000 site lattices. Extended cluster distribution tables are given in the supplementary material. (See Supplementary Material Available paragraph at the end of this article.) We see that as the size of the lattice increases we can get a better linear fit at concentrations closer to the critical point. We note also that the value of the slope is strongly dependent on the value of c,. Table I clearly illustrates this fact. For example, using 0.5000,0.5001, and 0.5002 as c, for a triangular 32000 X 32000 lattice, we get for y the values 2.457,2.404, and 2.351, respectively, implying that an 0.05 uncertainty in y yields an 0.0001 uncertainty in c,. We observe also from Table I and the monolayer data of Table I1 that y increases as the lattice size increases. Up to this point, we assumed that we know c, and calculated y. Now we take the opposite approach and calculate c, by using the theoretical value of y for two-dimensional lattices, y = 43/18. Using our simulation result and the theoretical y, we get that c, lies between 0.5000 and 0.5001. Using the theoretical y for the square lattice for a lattice of 32000 X 32000 sites, we see that by using y values from Table y out
X 400
sites,
TABLE 11: y for Square Lattices
edge length 400
lavers 1 2 4
8 16 32 400 4 000
16 000 32 000
1 2 3 4 1
1 1
1
c,
y
0.5927 0.475 0.398 0.3561 0.331 0.3218 0.3135 0.5927 0.4758 0.4271 0.4008 0.5931 0.5932 0.5931 0.5930
2.17 2.30 2.30 2.17 1.95 1.87 1.78 2.18 2.33 2.37 2.36 2.389 2.428 2.400 2.371
correlation configurations 0.99993 0.99991 0.99996 0.99986 0.99999 0.99995 0.99981 0.99990 0.99992 0.99998 0.99993 0.99989 0.999979 0.999990 0.999997
50 39 51
5 5 5 2 10 1
1 1 1 1 1 1
I1 lying between 2.400 and 2.371 we get 0.5930 Ic, I0.5931. However, using the "finite lattice" y from Table I (2.351), we get from Table I1 that c, = 0.5929 0.0002. It is rather clear that the Monte Carlo method can only yield semiquantitative results for the y exponent even for very large
*
The Journal of Physical Chemistry, Vol. 91, No. 1, 1987 221
Refined Monte Carlo Simulations of Static Percolation 2. 5
5
2. 3
1000
10
U (Y 0
rc
0
2
\\
2. 1
u 0
0
100
N
L
> H
1.9
10 d
1.7
1
100
10
1000
1
1
Figure 2. Effect of multiple layers of square lattices on y. y vs. the number of layers. (Edge length: 400 (-), 4000 (- -), 16000 (- -), and 32 000 (-- -).
400
1000 2000 4000
Wl
0.0055 0.0046 0.0034
0.0017 0.001 1
configurations 50 50 30 14 17
lattices. Obviously, when there are no alternate methods available for determining y, the Monte Carlo method is indispensable. In those cases the Monte Carlo can provide an estimate for the value of y. For example, we were able to investigate the change in y as we go from a single layer, two-dimensional rectangular 4000 X 4000 lattice to an n layer 400 X 400 X n cubic system. We see that y gradually decreases as the number of layers increases as illustrated in Table I1 and Figure 2. The apparent maximum in y for two layers for 400 X 400 in Figure 2 is an artifact of the small size of these lattice. For the larger lattices, 16000 X 16000 and 32 000 X 32 000, this “anomaly” disappears. We note that k i= 1 over a wide concentration range (0.59 to 0.01).
Correlation Length Exponent u To obtain the correlation length exponent u from the Monte Carlo data, we have used the Levinshtein et al.” method for u. Levinshtein et al. evaluate u from the dispersion of the values of c, vs. the linear dimension 1 of the simulated lattice. The Monte Carlo method appears to be a natural method for determining u since the outcome of multiple Monte Carlo runs provides a dispersion of random variables. The dispersion Wl for a lattice of linear dimension 1 is given by W/ = ( ( c , -
(C,)l2)
100
1000
0000
Figure 3. v for a square lattice. A log-log plot of inverse of dispersion vs. the lattice edge length.
W e assume that the greater scatter in our results is due to our smaller sample size.
TABLE III: Dispersion of c. with Lattice Size
edge length 256
10
Lottlcc Edge., 1
Number of Layers
(4)
where ( ) denotes a mean value for all simulated lattice configurations of linear dimension 1. W, is given as a function of 1 W1 -- Bl-(I/Y) (5) where B is a constant. Our results for W/vs. 1 for large lattices are listed in Table 111. Our five data points for large lattices plus the results of Levinshtein et al. for small lattices are also displayed in Figure 3. The value of u is determined by Levinshtein et al. to be 1.33 f 0.04 for 1 values between 4 and 128. When the W/values from Table I11 are used together with 1 = 64 and 128, we obtain u = 1.34 f 0.08. This is in agreement with the v obtained for the small lattices. (11) Levinshtein, M. E.; Shklovskii, B. I.; Shur, M. S.; Efros A. L. Sou. Phys. JETP 1976, 42, 197.
Discussion There is a range of literature vlaues of c, for the square lattice.’-*I The question we raise is as follows: for a finite lattice, is there a single value for c,, or will it depend on the criterion used for obtaining it (in addition to depending on the size, boundary conditions, and possibly the random number generator) ? Among the criteria used were boundary to boundary connectivities,” derivative of the second moment of the cluster distribution,12 and reduced average cluster size (see (2)).8-10 Nakanishi and StanleyI3 and Hoshen, Stauffer, et aI.l4 used the criterion y = 7’. Reynolds et al.15used a similar criterion of a c, that is symmetric with respect to position space renormalization group extrapolations. Quinn et a1.16 used the criterion of best fit of cluster groups to the scaling relation for 6. RapaportI7 used the latter criterion as well as a criterion based on integrated cluster distributions. Recently, a number of new criteria have been added.’8-20 The actual values for c, for the square lattice have oscillated for the last few yearss,14-20between 0.5927 and 0.5931. More interesting, probably, is c, for the “real” (finite) triangular lattice (as compared to the theoretical 0.5). Hoshen et al.I4 suggested 0.500 85, for a 4000 X 4000 lattice with free boundary conditions. Our results here suggest a value 0.5000 Ic, I0.5001 for our 32 000 X 32000 lattice which is consistent. By the same token, our value for the square lattice, c, = 0.5929 f 0.0002, is on the high side of recent “oscillations”. It is “below” the value of Reynolds et al.ls of 0.5931 f 0.0006 and “above” the value of 0.5927 f 0.0003 suggested by Hoshen et aL8 and in the most recent work of Rapaport,” 0.5927 f 0.0001, and the Stauffer Monograph7 quoted value of 0.59275. Another recent workI9 gives c, = 0.592802 f 0.00001 derived by a diffusion front model, and simulation results giving 0.592 745 were obtained very recently (12) Dean, P.; Bird, N. F. Public. No. Ma61, National Physics Laboratory Mathimatical Division, 1966. (13) Nakanishi, H.; Stanley, H. E. J. Phys. A: Math. Gen. 1978, 1 2 , L189. (14) Hoshen, J.; Stauffer, D.; Bishop, G . H.; Harrison, R. J.; Quinn, G. D.J . Phys. A . Math. Gem 1979, 12, 1285. (15) Reynolds, P. J.; Stanley, H. E.; Klein, W. Phys. Rev. B 1980, 21, 1223. (16) Quinn, G. D.; Bishop, G. H.; Harrison R.J. J . Phys. A: Math. Gen. 1976, 9, L9. (17) Rapaport, D. C. J. Phys. A 1985, 18, L175. (18) Gebele, T. J . Phys. A 1984, 17, L51. (19) Rosso, M.; Gouyet, F.; Sapoval, B. Phys. Rev. B 1985, 32, 6053. (20) Ziff, R. M. Phys. Reo. Lett. 1986, 56, 545, and private communica-
tion.
(21) Newhouse, J. S. Ph.D. Thesis, The University of Michigan, 1985.
222
J. Phys. Chem. 1987, 91, 222-226
by Ziff,” based on cluster perimeters. The question arises whether the low error bars claimed in recent references can be trusted. If yes, then the real question that arises is how close an agreement can we expect, for finite lattices, from different criteria? Our v value is consistent with previous values’ as well as the theoretical conjecture, v = 4/3. Within our error limits, the scaling hypothesis (combined with the den Nijs conjecture) is consistent for our finite, free-boundary lattice. Our layer problem studies provided few surprises. y is monotonic with the number of layers. However, we emphasize that we preserved, roughly, the total number of sites. (However, we should note that a single 400 X 400 square lattice simulation gave2* lower y values, leading to an apparent nonmonotonic relationship between y and the number of layers.) The value of y = 1.78 for the cubic lattice (Table 11, 400 X 400 X 400) is consistent with the literature’ value of 1.8 (however, our c, = 0.3135 differs from the literature value of 0.31 17). Why the preexponent k should be unity over a wide concentration range is not clear to us.
Conclusion (1) y is very sensitive to c,. (2) Knowledge of y improves c,. (3) The method (HK) of c, based on the maximum in the mean cluster size (excluding the largest cluster) is consistent with the determination of c, using the theoretical value of y. (4) Simulations on a lattice of IO9 sites gives c, values to an accuracy of about using theoretical y values. Square and triangular lattices of the size lo9 give y values consistent with universality (independence of y on topology) and with the theoretical values to within 2%. (5) Lattices up to lo9 sites give v values within 1% of the theoretical values and the universality hypothesis. (6) Effective y exponents for multiple layers show a monotonic crossover from the two-dimensional to the three-dimensional values.
Supplementary Material Available: Extended cluster distribution tables (69 pages). Ordering information is given on any current masthead page.
Behavior of Ethanol in Varlous Binary Solutions: Dlfference Raman Spectroscopy on the C-H Stretching Vibrations K. Kamogawa, S. Kaminaka, and T. Kitagawa* Institute f o r Molecular Science, Okazaki National Research Institutes, Myodaiji, Okazaki, 444 Japan (Received: June 30, 1986)
Frequency shifts of the C-H stretching modes of ethanol upon mixing with various liquids were observed as a function of interaction the concentration,and their shift behaviors were analyzed by using the homogemus (Avdand heterogeneous (bm) factors proposed previously. The values of AvABat x = 1.0, where x is the mole fraction of ethanol, became smaller as the polarizability of the perturbing species became larger. This suggested that the dispersion force among various intermolecular forces plays an essential role in determining the frequency shifts of the C-H stretching modes upon mixing. The AvAA vs. x curve of the C2H,0H/H20 solution exhibited noticeable similarity to the correspondingcurve of the partial molar volume of ethanol, as anticipated from the previous study on the CH3CN/H20and (CH3)2CO/H20solutions. This study provides spectroscopic evidence for the presence of discrete composition segments in the binary solutions.
An aqueous solution of the organic molecule having both hydrophobic and hydrophilic groups seems to consist of some kind of clusters, the number of which is a smooth function of the concentration.’S2 However, there is a possibility that a change of the concentration does not mean a simple population change of the components but involves qualitative alteration of the solution structure. This idea was recently incorporated in t e m of “discrete composition segment^"*^ to explain the excess molar volume upon mixing. It is desirable to spectroscopically demonstrate the presence of the discrete composition segments, that is, a solution structure specific to individual ranges of the concentration. For this purpose Raman spectroscopy seems promising; since molecular vibrations are generally sensitive to intermolecular interactions, vibrational frequencies are expected to explore a change of the microscopic environments around a probe molecule in a solution. We have constructed a difference Raman spectrometer6 and in fact observed the concentration-dependent frequency shifts of the C-H stretching modes upon mixing of two l i q ~ i d s . ~ ~ ~ (1) Marcus, Y. In Introduction to Liquid State Physics; Wiley: London, 1977; Chapter 5. (2) Larkin, J. A. J. Chem. Thermodyn. 1975, 7, 137. (3) Davis, M. I. Thermochim. Acta 1983, 63, 67. (4) Davis, M. I. Thermochim. Acta 1983, 71, 59. (5) Davis, M. I. Thermochim. Acta 1984, 81, 137. (6) Kamogawa, K.; Tajima, K.; Hayakawa, K.; Kitagawa, T. J . Phys. Chem. 1984, 88, 2494. (7) Kamogawa, K.; Kitagawa, T. J . Phys. Chem. 1985, 89, 1531
0022-3654/87/2091-0222%01.50/0
In the analysis of the frequency shifts for the binary solutions, we noticed that the homogeneous and heterogeneous interaction factors derived from the observed frequency shifts scale the magnitudes of the intermolecular interactions between molecules of the same species and those between molecules of different species, respectively, and that each quantity reflects a microscopic association ~ t r u c t u r e . ~These , ~ quantities have been obtained for the aqueous solutions of methanol, acetonitrile, and acetone, and all of them suggested the transition from the microscopic liquid phase to a cluster phase upon dilution by water. It was also stressed that methanol behaved differently from acetone and acetonitrile in a low-concentration region. The peculiarity might be ascribed to insufficient hydrophobicity of the methyl group of methanol. To examine this difference as well as an empirical rule about the shift factors and thermodynamic properties proposed previously, ethanol as a more hydrophobic alcohol9 was investigated in the present study. Experimental Section The Raman difference spectrometer consists of a single monochromator (Jasco, CT-50), an intensified diode array detector (PAR 1420), and a microcomputer (Oki, if800). Instrumental details were described elsewhere.6 The 514.5-nm line (200 mW) (8) Kamogawa, K.; Kitagawa, T. J. Phys. Chem. 1986, 90, 1077. (9) Rowlinson, J. S. In Liquid and Liquid Mixtures, 2nd ed.;Butterworth: London, 1977; p 168.
0 1987 American Chemical Society
