- Email: [email protected]
A fast cluster counting algorithm for percolation on and off lattices
Computer Physics Communications ELSEVIER
Computer PhysicsCommunications 109 (1998) I-5
A fast cluster counting algorithm for percolation on and off lattices E. S t o l l i Physik-lnstimt, Universitiit Ziirich-lrchel. Winterthure~tr. 190, CH-8057 Ziirich. Switzerland
Received 5 November 1997; revi.~l I December 1997
Abstract A fast counting algorithm is used to enumerate the clusters for determining their statistics, for investigation the percolation behavior of the largest cluster and for calculating the fractai dimensions D. The algorithm is suitable for parallel processing. It can be expanded for off-lattice systems, where particles are considered to be connected when they are separated by less than a given threshold. For determining the fractal dimensions D, the procedure is very fast in investigating supercells, where replica of smaller boxes with off-lattice structures are arranged. © 1998 Elsevier Science B.V. PACS: 02.50.Ey; 02.7u.Lq; 47.53.+n; 61.43.Hv; 61.20Ja; 64.60.Ak; 64.60.-i Keywords: Stochasticprocesses;Monte Carlo: Fractals;Computer simulation; Percolation;Phase transitions; Droplet model
1. Introduction In the past, a variety of studies have been performed to investigate cluster statistics [ I - 3 ] , percolation thresholds [2-6] and fractal properties [4-6]. A fast cluster counting algorithm was established to determine the critical percolation concentration for the site percolation problem on lattices [7]. To generalize the percolation and fractal problem for off-lattice systems an alternative method for cluster counting is needed. Similar basic ideas for enumerating the clusters in lattices were developed at the same time as in Ref. [7]. In a finite system with periodic boundary conditions, percolation occurs when the largest cluster is connected across the whole system. This connection property is investigated. Such a direct test is unavoidable when such a system is used to
calculate the dynanlics of fractal systems [6]. To determine the dynamics of the excitations in fractal systems, the fractons, dynamic matrices are generated and their eigenvalues and eigenvectors are calculated [6]. For large systems the computer code can be parallelized by dividing the whole space into smaller subregions. Such algorithms using the method of Ref. [7] are described in Refs. [8-10]. These algorithms are improved by the authors of Refs. [ 11 ] and [ 12] who presented very sophisticated procedures to scatter and gather the cluster information over and from a large amount of distributed memory processors. With our counting algorithm this parallelization can be simplified. When elements not arranged on a lattice belong to a cluster, they are defined as connected when they are separated by less than a given threshold length [ 13]. From this ensemble of bonds and their accompanying sites, subsets can be selected by random number gen-
l Correspondingauthor: E-mail:[email protected]. 0010-4655/98/$19.00 ~) 1998 Elsevier Science B.V. All rightsreserved. PII S0010-4655 (97) 00142-2
E. Stoll/Computer Physics Communications 109 (1998) I-5 erators for bond or site percolation. Percolation clusters and fractal exponents have also been determined in off-lattice fractal systems [ 13,14]. This selective choice of bonds or sites protects the system against multifractality [ 15 ].
2. Theory and informatics requirements In the cluster counting procedure each element and its environment are investigated for connections. All connected elements belong m clusters. All clusters are enumerated. When two connected elements are not enumerated, the largest order number I of all clusters is increased by one to l + I and given to both elements. When one element has not been enumerated, this element takes the order number of the other one. When both elements of this cluster have the same order number, no action is required. Things become more difficult when both elements belong to two clusters with different order numbers i and j , and i < j. To solve this problem, index registers are used, as was oone in Ref. [7]. Prior to applying the cluster counting procedure, all elements in a first register are put to zero, and all elements Mk in a second register are assigned with their order number k. In the first register, the ith element N/becomes j . In the second register, the jth element M j becomes i. This arrangement can easily be expanded for a concatenation of more clusters, NJ, = J2, Ivi2 = J3 , " " , Nj .... = j,,, N j,, = 0, forjl
(I)
and M,, = M j .... = . . . .
Mj, = j ~ .
(2)
When the order numbers i ~ j but Mi = M i, both clusters are already concatenated, and no further action is required. Elsewhere, for Mi = il, M j = Jl, and il < j r , both registers have to be updated. Because both chains {Ni} and {Nj} are monotonically arranged (see Eq. ( I ) ) , the merging of these chains and the updating of all elements of the chain {Mj} with il can easily be performed: Nij contains i2. As long as, for m > 2, all i,n are positive and i,n < Jl, Ni.._, remains i,n. When i m > j l , Ni._, becomes jr, and Mj, becomes il. Afterwards, as long as, for n > 2,
all positive j , < i,,, Nj,,_~ remains j,, and Mj,, becomes il. Elsewhere, Nj,,_, becomes in,, and this "'changing procedure" is continued until either j , becomes zero or i,, becomes zero. When jn becomes zero first then Nj,,_, becomes in,, Mj,,_, becomes il, and the updating process is completed. When, on the other hand, in, becomes zero first, Ni.,_, becomes jn. Afterwards, the remaining chain of Mg,, becomes i~ until jn becomes zero. The contents of the first and second registers now satisfy the requirements of Eqs. ( I ) and (2). Unfortunately, for use with today's last vector computers, the recursive nature of this updating algorithm does not allow computer code vectorization. After completing the cluster enumeration and concatenation procedure for all elements in a given region, the content of the second register is modified. Mi= j = l = I. j is increased monotonically. When Mj = j , I is increased by one and M j takes the value/. Elsewhere, the value Mj is changed to M~,. After performing this modification, ! corresponds to the number of all clusters in the system with more than one element. The elements of all those clusters are renumbered by replacing their order numbers j by Mj. To be ready for further cluster numbering and concatenating steps, all elements in the first register are put to zero, and all elements My in the second register are again assigned with their order number j. For very large systems, the concatenated chains in both registers become rather long during the counting period. Because the updating algorithm, called by every cluster amalgamation, must run in both registers through one of these chains, the computing speed slows down. Therefore, for very large and multi-concatenated systems it is more favorable to subdivide the whole system into smaller subsystems and to renumber all elements of all considered systems prior applying the counting procedure on the next subsystem. The next step to save time is to parallelize the procedure and to distribute the computing load on several processors. This was proposed and implemented in Refs. [ 8-12] without using the information of our second register. The system is again subdivided into smaller subregions. In each subregion k the cluster counting procedure is applied for all elements within this region. The elements then have to be renumbered. To avoid ambiguities, the order numbers of the clusters in neighboring subregions must be different. The
E. Stofl/Computer Physics Communications 109 (1998) I-5
contents of all M~ have to be modified to MJ' with k-I .!' = M jk + E l i M~: i=1
(3)
where li is the number of clusters in the ith subregion. Afterwards, the cluster counting procedure will be continued over all common border area, lines and corners of all subsystems as it was reported in Refs. [ ! I ] and [ 12 ], and the elements are again renumbered. The application of the information of the second registers in each processor to compare the cluster order numbers can help to reduce the communication complexity in the already published procedures [ 1 !,12]. For simulation of random fractals [5] and their dynamics [6] it is indispensable to know whether or not the largest cluster is connected across the whole system. This problem can be solved by the direct comparison of the cluster order numbers of the elements on the opposite system bord =s when this cluster ?ges not wind like a snake through the periodic borders, as shown in Fig. I. This cluster seems to be disconnected using the oversimplified arguments mentioned above. To avoid this difficulty, further effort is needed. All particles only connected via periodic boundaries are marked with order numbers f and ! + !. Across the borders not considered for percolation tests, the concatenation is marked in the first and second register, according to Eqs. ( I ) and (2). W~'en the percolating cluster is connected directly via the periodic border, particles opposite the border belonging to this cluster have an identical cluster order number M~ and Mj. This number is noted in the percolation cluster register. When Mi --/: M j, Mi and M) are saved in a third and fourth register. After the investigation of all sites on both sides of the periodic border, the contents of the third and fourth registers are compared and the connection is marked, in well prepared fifth and sixth registers, as Mta, and Mtaj+t according to Eqs. ( 1) and (2). (l is the number of all clusters with more than one element in the system.) MMj+t is used instead of MM, to avoid trivial identities for neighboring elements opposite to the periodic border. When, for an element Mi in the third register, the values MM, and M~,+t are identical in the sixth register, the cluster with the order number Mi is connected, and Mi is also saved in the percolation cluster register.
,~ o°ggggo o°g~ggo hg. . . . . . . . gggoggggggg° mm.g. . . . . . . goo o~ o
oogggg°
:
i °°ggggo
gg=:::::;
o.
Fig, I. A percolation cluster winding like a snake tl'a'onghthe horizontal cluster border. The different symbols denote different clusterorder numbersas they are generatedpriorto the percolation cluster test. Using periodic boundary conditions, the cluster distribution and the existence of percolation clusters are translation invariant. Our procedures have been tested using this invariance: a large number of systems with identical distributions of occupied sites or bonds were generated and then shifted over an arbitrary number of sites in each space dimension. As identical cluster distributions and identical percolation properties of the largest clusters were achieved for identical original site or bond distributions, we may then conclude that our programs passed this test successfully.
3. lsotropic
systems
(continuum
percolation)
lsotropic systems can be generated off-lattice with isotropically iateracting particles enclosed in boxes [ 13,14]. The long range order of such fluids breaks down and only a short range order exists. Its two point correlation shows that at a certain distance nearly all nearest neighbor connections are included and nearly all other connections are excluded. We defined a threshold length 2.5% larger than the potential minimum between two particles. Particles that are separated less than this threshold are assumed as connected. A careful investigation of such Lennard-Jones systems showed that the n u m ~ r of bonds using this procedure is rather insensitive to small fluctuations of this threshold value [ 13,14]. Each of those thus
E. Stoll/Computer Physics Communications 109 (1998) I-5
...... /
.....
~!:: .........
:~
l--I~"o o
° o-.-2~
'I ........ / : .........
......
::.
~::i.......... i i
Fig. 2. Connection scheme in two dimensions. The arrows denote the three different types: I. both particles are in the bulk, 2. the
particles are connected via the area border, 3. the connection is performed via the comer. For throe dimensionsthe scheme has to be expended accordingly determined bonds is stored together with the labels of both of its adjoining particles. The whole set of bonds is subdivided into different classes according to the connection paths. Some bonds connect both particles inside the simulated volume. Periodic boundary conditions are applied. The remaining bonds connect two particles only by crossing box borders. Concatenating the system with its replica, the bonds connect particles via the bounding surfaces, edges and corners of the original system as shown for two dimensions in Fig. 2. With the stored bcnds and the labels of both their adjoining particles, the cluster counting procedure of Section 2 is applied. Random numbers between zero and one are generated and compared with a supplied value p. When the random number is smaller than or equal to p, the considered bond exists, or the investigated site is occupied. In the cluster counting procedures only existing bonds are used for bond percolation. For site percolation the applied bonds have to connect two occupied sites. The different types of connections have also to be considered in the subroutines that investigate the percolation property of the largest cluster described in Section 2. The connections via the opposite surfaces and the joined edges and corners are introduced to include all possible configurations. The box-counting method has been applied for
boxes of different sizes [ 13,14], for the investigation of the fractai behavior, and for the determination of the fractal dimension D. In principle this can be performed by the simulation of a large number of independent spatial particle distributions and different random number chains for bond or site choices. Unfortunately, the molecular dynamical simulations of different particle distributions need the largest amount of the computer effort. Furthermore, the fluctuations are very large. Therefore, we simulated only a few different spatial distributions in systems which are not greatly expanded, and conne:ted an array of those boxes via the periodic boundary conditions in all space dimensions. The subdivision of the bonds into different classes according to the connection paths in the inner of the bulk or via the bounding surfaces, edges and corners in the concatenated system greatly simplifies the determination of the Nb bonds and both their accompanying sites. The number of operations is reduced from the order N~, to the order Nh. With the same spatial structure but with different random chains for the bond or site choices a lot of systems were simulated. The fractal dimension D could be determined [13,14]. Furthermore, it was possible to show that our off-lattice systems are fractal [ 13,14 I and multifractality [ 15] is not needed.
4. S u m m a r y and conclusions With our algorithm we are able to enumerate the clusters and to calculate their statistics. In contradiction to Ref. [7], two registers are used to readjust the numbers of concatenated clusters. Our first register corresponds to that used by Ref. [7], with the small exception that we have applied positive values ij+l tbr the contents of Nij in the inside of the chain and zero for the top chain element instead of the number of elements of the concatenated cluster. Unlbrtunately, with our algorithm we use slightly more computer time than with those of Ref. [7] when two clusters amalgamate, because we then have to update the chains in two registers. However, even by comparing elements of clusters that have been marked as concatenated, the easy comparison of their common cluster number via the second register can save a lot of computer eflbrt. Furthermore, our algorithm allows one to directly determine percolation clusters, to calculate the fractal
E. Stoll/Computer Physics Communications 109 (1998) 1-5
dimensions D, and by applying the second register to improve the parallel processing procedures described in Refs. [ ! I ] and [ 12]. It can be expanded to offlattice systems, where those particles are considered as connected when they are separated by less than a given threshold. For determining the fractal dimensions D, supercells with periodic continued off-lattice structures are used, and the computer effort increases only linearly as a function of the number of bonds or of considered particles.
Acknowledgements We would like to thank K. Binder, E. Courtens, C. Domb, Ch. Hanser, H. Mtiller-Krumbhaar, T. Schneider, J. Singer, D. Stauffer, and Ch. Stern for fruitful and constructi~,e discussions, and William Harris for proofreading the manuscript. Parts of the work were performed at the IBM Research Division, Zurich Research Laboratory in Riischlikon prior to the author's joining the University of Zurich.
References [ I ] E. Stoll, K. Binder, T. Schneider; Phys. Rev. B 6 (1972) 2777. [21 H. Miiller-Kmmbhaar, E.P. Stoll, J. Chem. Phys. B 65 (1976) 4294. 131 E. Stoll, C. Domb, J. Phys. A 12 (1979) 1843. 141 E.P. Stoll, M. Kolb, Physica A 185 (1992) 222. [5l E. Counens, R. Vacher, E. Stoll, Physica D 38 (1989) 41. 161 E. Stoll, M. Kolb. E. Courtens. Phys. Rev. Lea. 68 (1992,) 2472. [71 J. Hoshen. R. Kopelman. Phys. Rev. B 14 (1976) 3438. 181 A.N. Burkitt, D.W. Heermann, Compm. Phys. Commen. 54 (1989) 201. 191 D.W. Heennan,, A.N. Burkitt, Parallel Algorithms in Computational Science (Springer, Heidelberg, 1991). 110l R.C. Brower, P. Tamayo, B. York, J. Star. Phys. 63 (1991) 73. [111 M Bauemfeind. R. Hackl. H.-G. Matu, is, J. Singer, Th. Husslein. I. Morgenstem. Physica A 212 (1994) 277. [ 12] M. Flanigan P. Tamayo. Physica A 215 (1995) 461. [ 13] E. SIolL Ch. Stem, P. Stucki. Physica A 230 (1996) Ii. {141 E. Stoll, Ch. Stem, P. Stucki, Heir. Phys. Acta 69 (1996) 21. l l5l J.H.J. van Opheusden, M.T.A. Bos, G. van der Kaaden. Physica A 227 (1996) 183.
Computer PhysicsCommunications 109 (1998) I-5
A fast cluster counting algorithm for percolation on and off lattices E. S t o l l i Physik-lnstimt, Universitiit Ziirich-lrchel. Winterthure~tr. 190, CH-8057 Ziirich. Switzerland
Received 5 November 1997; revi.~l I December 1997
Abstract A fast counting algorithm is used to enumerate the clusters for determining their statistics, for investigation the percolation behavior of the largest cluster and for calculating the fractai dimensions D. The algorithm is suitable for parallel processing. It can be expanded for off-lattice systems, where particles are considered to be connected when they are separated by less than a given threshold. For determining the fractal dimensions D, the procedure is very fast in investigating supercells, where replica of smaller boxes with off-lattice structures are arranged. © 1998 Elsevier Science B.V. PACS: 02.50.Ey; 02.7u.Lq; 47.53.+n; 61.43.Hv; 61.20Ja; 64.60.Ak; 64.60.-i Keywords: Stochasticprocesses;Monte Carlo: Fractals;Computer simulation; Percolation;Phase transitions; Droplet model
1. Introduction In the past, a variety of studies have been performed to investigate cluster statistics [ I - 3 ] , percolation thresholds [2-6] and fractal properties [4-6]. A fast cluster counting algorithm was established to determine the critical percolation concentration for the site percolation problem on lattices [7]. To generalize the percolation and fractal problem for off-lattice systems an alternative method for cluster counting is needed. Similar basic ideas for enumerating the clusters in lattices were developed at the same time as in Ref. [7]. In a finite system with periodic boundary conditions, percolation occurs when the largest cluster is connected across the whole system. This connection property is investigated. Such a direct test is unavoidable when such a system is used to
calculate the dynanlics of fractal systems [6]. To determine the dynamics of the excitations in fractal systems, the fractons, dynamic matrices are generated and their eigenvalues and eigenvectors are calculated [6]. For large systems the computer code can be parallelized by dividing the whole space into smaller subregions. Such algorithms using the method of Ref. [7] are described in Refs. [8-10]. These algorithms are improved by the authors of Refs. [ 11 ] and [ 12] who presented very sophisticated procedures to scatter and gather the cluster information over and from a large amount of distributed memory processors. With our counting algorithm this parallelization can be simplified. When elements not arranged on a lattice belong to a cluster, they are defined as connected when they are separated by less than a given threshold length [ 13]. From this ensemble of bonds and their accompanying sites, subsets can be selected by random number gen-
l Correspondingauthor: E-mail:[email protected]. 0010-4655/98/$19.00 ~) 1998 Elsevier Science B.V. All rightsreserved. PII S0010-4655 (97) 00142-2
E. Stoll/Computer Physics Communications 109 (1998) I-5 erators for bond or site percolation. Percolation clusters and fractal exponents have also been determined in off-lattice fractal systems [ 13,14]. This selective choice of bonds or sites protects the system against multifractality [ 15 ].
2. Theory and informatics requirements In the cluster counting procedure each element and its environment are investigated for connections. All connected elements belong m clusters. All clusters are enumerated. When two connected elements are not enumerated, the largest order number I of all clusters is increased by one to l + I and given to both elements. When one element has not been enumerated, this element takes the order number of the other one. When both elements of this cluster have the same order number, no action is required. Things become more difficult when both elements belong to two clusters with different order numbers i and j , and i < j. To solve this problem, index registers are used, as was oone in Ref. [7]. Prior to applying the cluster counting procedure, all elements in a first register are put to zero, and all elements Mk in a second register are assigned with their order number k. In the first register, the ith element N/becomes j . In the second register, the jth element M j becomes i. This arrangement can easily be expanded for a concatenation of more clusters, NJ, = J2, Ivi2 = J3 , " " , Nj .... = j,,, N j,, = 0, forjl
(I)
and M,, = M j .... = . . . .
Mj, = j ~ .
(2)
When the order numbers i ~ j but Mi = M i, both clusters are already concatenated, and no further action is required. Elsewhere, for Mi = il, M j = Jl, and il < j r , both registers have to be updated. Because both chains {Ni} and {Nj} are monotonically arranged (see Eq. ( I ) ) , the merging of these chains and the updating of all elements of the chain {Mj} with il can easily be performed: Nij contains i2. As long as, for m > 2, all i,n are positive and i,n < Jl, Ni.._, remains i,n. When i m > j l , Ni._, becomes jr, and Mj, becomes il. Afterwards, as long as, for n > 2,
all positive j , < i,,, Nj,,_~ remains j,, and Mj,, becomes il. Elsewhere, Nj,,_, becomes in,, and this "'changing procedure" is continued until either j , becomes zero or i,, becomes zero. When jn becomes zero first then Nj,,_, becomes in,, Mj,,_, becomes il, and the updating process is completed. When, on the other hand, in, becomes zero first, Ni.,_, becomes jn. Afterwards, the remaining chain of Mg,, becomes i~ until jn becomes zero. The contents of the first and second registers now satisfy the requirements of Eqs. ( I ) and (2). Unfortunately, for use with today's last vector computers, the recursive nature of this updating algorithm does not allow computer code vectorization. After completing the cluster enumeration and concatenation procedure for all elements in a given region, the content of the second register is modified. Mi= j = l = I. j is increased monotonically. When Mj = j , I is increased by one and M j takes the value/. Elsewhere, the value Mj is changed to M~,. After performing this modification, ! corresponds to the number of all clusters in the system with more than one element. The elements of all those clusters are renumbered by replacing their order numbers j by Mj. To be ready for further cluster numbering and concatenating steps, all elements in the first register are put to zero, and all elements My in the second register are again assigned with their order number j. For very large systems, the concatenated chains in both registers become rather long during the counting period. Because the updating algorithm, called by every cluster amalgamation, must run in both registers through one of these chains, the computing speed slows down. Therefore, for very large and multi-concatenated systems it is more favorable to subdivide the whole system into smaller subsystems and to renumber all elements of all considered systems prior applying the counting procedure on the next subsystem. The next step to save time is to parallelize the procedure and to distribute the computing load on several processors. This was proposed and implemented in Refs. [ 8-12] without using the information of our second register. The system is again subdivided into smaller subregions. In each subregion k the cluster counting procedure is applied for all elements within this region. The elements then have to be renumbered. To avoid ambiguities, the order numbers of the clusters in neighboring subregions must be different. The
E. Stofl/Computer Physics Communications 109 (1998) I-5
contents of all M~ have to be modified to MJ' with k-I .!' = M jk + E l i M~: i=1
(3)
where li is the number of clusters in the ith subregion. Afterwards, the cluster counting procedure will be continued over all common border area, lines and corners of all subsystems as it was reported in Refs. [ ! I ] and [ 12 ], and the elements are again renumbered. The application of the information of the second registers in each processor to compare the cluster order numbers can help to reduce the communication complexity in the already published procedures [ 1 !,12]. For simulation of random fractals [5] and their dynamics [6] it is indispensable to know whether or not the largest cluster is connected across the whole system. This problem can be solved by the direct comparison of the cluster order numbers of the elements on the opposite system bord =s when this cluster ?ges not wind like a snake through the periodic borders, as shown in Fig. I. This cluster seems to be disconnected using the oversimplified arguments mentioned above. To avoid this difficulty, further effort is needed. All particles only connected via periodic boundaries are marked with order numbers f and ! + !. Across the borders not considered for percolation tests, the concatenation is marked in the first and second register, according to Eqs. ( I ) and (2). W~'en the percolating cluster is connected directly via the periodic border, particles opposite the border belonging to this cluster have an identical cluster order number M~ and Mj. This number is noted in the percolation cluster register. When Mi --/: M j, Mi and M) are saved in a third and fourth register. After the investigation of all sites on both sides of the periodic border, the contents of the third and fourth registers are compared and the connection is marked, in well prepared fifth and sixth registers, as Mta, and Mtaj+t according to Eqs. ( 1) and (2). (l is the number of all clusters with more than one element in the system.) MMj+t is used instead of MM, to avoid trivial identities for neighboring elements opposite to the periodic border. When, for an element Mi in the third register, the values MM, and M~,+t are identical in the sixth register, the cluster with the order number Mi is connected, and Mi is also saved in the percolation cluster register.
,~ o°ggggo o°g~ggo hg. . . . . . . . gggoggggggg° mm.g. . . . . . . goo o~ o
oogggg°
:
i °°ggggo
gg=:::::;
o.
Fig, I. A percolation cluster winding like a snake tl'a'onghthe horizontal cluster border. The different symbols denote different clusterorder numbersas they are generatedpriorto the percolation cluster test. Using periodic boundary conditions, the cluster distribution and the existence of percolation clusters are translation invariant. Our procedures have been tested using this invariance: a large number of systems with identical distributions of occupied sites or bonds were generated and then shifted over an arbitrary number of sites in each space dimension. As identical cluster distributions and identical percolation properties of the largest clusters were achieved for identical original site or bond distributions, we may then conclude that our programs passed this test successfully.
3. lsotropic
systems
(continuum
percolation)
lsotropic systems can be generated off-lattice with isotropically iateracting particles enclosed in boxes [ 13,14]. The long range order of such fluids breaks down and only a short range order exists. Its two point correlation shows that at a certain distance nearly all nearest neighbor connections are included and nearly all other connections are excluded. We defined a threshold length 2.5% larger than the potential minimum between two particles. Particles that are separated less than this threshold are assumed as connected. A careful investigation of such Lennard-Jones systems showed that the n u m ~ r of bonds using this procedure is rather insensitive to small fluctuations of this threshold value [ 13,14]. Each of those thus
E. Stoll/Computer Physics Communications 109 (1998) I-5
...... /
.....
~!:: .........
:~
l--I~"o o
° o-.-2~
'I ........ / : .........
......
::.
~::i.......... i i
Fig. 2. Connection scheme in two dimensions. The arrows denote the three different types: I. both particles are in the bulk, 2. the
particles are connected via the area border, 3. the connection is performed via the comer. For throe dimensionsthe scheme has to be expended accordingly determined bonds is stored together with the labels of both of its adjoining particles. The whole set of bonds is subdivided into different classes according to the connection paths. Some bonds connect both particles inside the simulated volume. Periodic boundary conditions are applied. The remaining bonds connect two particles only by crossing box borders. Concatenating the system with its replica, the bonds connect particles via the bounding surfaces, edges and corners of the original system as shown for two dimensions in Fig. 2. With the stored bcnds and the labels of both their adjoining particles, the cluster counting procedure of Section 2 is applied. Random numbers between zero and one are generated and compared with a supplied value p. When the random number is smaller than or equal to p, the considered bond exists, or the investigated site is occupied. In the cluster counting procedures only existing bonds are used for bond percolation. For site percolation the applied bonds have to connect two occupied sites. The different types of connections have also to be considered in the subroutines that investigate the percolation property of the largest cluster described in Section 2. The connections via the opposite surfaces and the joined edges and corners are introduced to include all possible configurations. The box-counting method has been applied for
boxes of different sizes [ 13,14], for the investigation of the fractai behavior, and for the determination of the fractal dimension D. In principle this can be performed by the simulation of a large number of independent spatial particle distributions and different random number chains for bond or site choices. Unfortunately, the molecular dynamical simulations of different particle distributions need the largest amount of the computer effort. Furthermore, the fluctuations are very large. Therefore, we simulated only a few different spatial distributions in systems which are not greatly expanded, and conne:ted an array of those boxes via the periodic boundary conditions in all space dimensions. The subdivision of the bonds into different classes according to the connection paths in the inner of the bulk or via the bounding surfaces, edges and corners in the concatenated system greatly simplifies the determination of the Nb bonds and both their accompanying sites. The number of operations is reduced from the order N~, to the order Nh. With the same spatial structure but with different random chains for the bond or site choices a lot of systems were simulated. The fractal dimension D could be determined [13,14]. Furthermore, it was possible to show that our off-lattice systems are fractal [ 13,14 I and multifractality [ 15] is not needed.
4. S u m m a r y and conclusions With our algorithm we are able to enumerate the clusters and to calculate their statistics. In contradiction to Ref. [7], two registers are used to readjust the numbers of concatenated clusters. Our first register corresponds to that used by Ref. [7], with the small exception that we have applied positive values ij+l tbr the contents of Nij in the inside of the chain and zero for the top chain element instead of the number of elements of the concatenated cluster. Unlbrtunately, with our algorithm we use slightly more computer time than with those of Ref. [7] when two clusters amalgamate, because we then have to update the chains in two registers. However, even by comparing elements of clusters that have been marked as concatenated, the easy comparison of their common cluster number via the second register can save a lot of computer eflbrt. Furthermore, our algorithm allows one to directly determine percolation clusters, to calculate the fractal
E. Stoll/Computer Physics Communications 109 (1998) 1-5
dimensions D, and by applying the second register to improve the parallel processing procedures described in Refs. [ ! I ] and [ 12]. It can be expanded to offlattice systems, where those particles are considered as connected when they are separated by less than a given threshold. For determining the fractal dimensions D, supercells with periodic continued off-lattice structures are used, and the computer effort increases only linearly as a function of the number of bonds or of considered particles.
Acknowledgements We would like to thank K. Binder, E. Courtens, C. Domb, Ch. Hanser, H. Mtiller-Krumbhaar, T. Schneider, J. Singer, D. Stauffer, and Ch. Stern for fruitful and constructi~,e discussions, and William Harris for proofreading the manuscript. Parts of the work were performed at the IBM Research Division, Zurich Research Laboratory in Riischlikon prior to the author's joining the University of Zurich.
References [ I ] E. Stoll, K. Binder, T. Schneider; Phys. Rev. B 6 (1972) 2777. [21 H. Miiller-Kmmbhaar, E.P. Stoll, J. Chem. Phys. B 65 (1976) 4294. 131 E. Stoll, C. Domb, J. Phys. A 12 (1979) 1843. 141 E.P. Stoll, M. Kolb, Physica A 185 (1992) 222. [5l E. Counens, R. Vacher, E. Stoll, Physica D 38 (1989) 41. 161 E. Stoll, M. Kolb. E. Courtens. Phys. Rev. Lea. 68 (1992,) 2472. [71 J. Hoshen. R. Kopelman. Phys. Rev. B 14 (1976) 3438. 181 A.N. Burkitt, D.W. Heermann, Compm. Phys. Commen. 54 (1989) 201. 191 D.W. Heennan,, A.N. Burkitt, Parallel Algorithms in Computational Science (Springer, Heidelberg, 1991). 110l R.C. Brower, P. Tamayo, B. York, J. Star. Phys. 63 (1991) 73. [111 M Bauemfeind. R. Hackl. H.-G. Matu, is, J. Singer, Th. Husslein. I. Morgenstem. Physica A 212 (1994) 277. [ 12] M. Flanigan P. Tamayo. Physica A 215 (1995) 461. [ 13] E. SIolL Ch. Stem, P. Stucki. Physica A 230 (1996) Ii. {141 E. Stoll, Ch. Stem, P. Stucki, Heir. Phys. Acta 69 (1996) 21. l l5l J.H.J. van Opheusden, M.T.A. Bos, G. van der Kaaden. Physica A 227 (1996) 183.