# An improved Rosenbluth Monte Carlo scheme for cluster counting and lattice animal enumeration

## 1 Introduction

The enumeration of lattice animals is an important problem in a variety of physical problems including nucleation , percolation  and branched polymers . A lattice animal is a cluster of $`N`$ connected sites on a lattice with given symmetry and dimensionality and we seek to enumerate all distinct animals with a given number of sites. Exact enumeration has been carried out for small lattice animals using a variety of methods  but the methods become computationally prohibitive for large animals. Many techniques have been used to enumerate larger lattice animals including various Monte Carlo growth schemes , a constant fugacity Monte Carlo method , an incomplete enumeration method  and reaction limited cluster-cluster aggregation .

In the following paper we describe an improvement of a method proposed by one of the authors  which was based on an extension of the scheme proposed by Rosenbluth and Rosenbluth  for enumerating self avoiding polymer chains. The central problem in using the Rosenbluth scheme for lattice animal enumeration is calculating the degeneracy of the clusters which are generated. In the method proposed by Care, the cluster growth was modified in a way which forced the degeneracy to be $`N!`$ where $`N`$ is the number of sites occupied by the lattice animal. However the resulting algorithm was fairly complicated to implement. An alternative method of correcting for the degeneracy had been proposed by Pratt . In this latter scheme the correcting weight is more complicated to determine and must be recalculated at each stage of the cluster growth if results are sought at each cluster size. However the Pratt scheme does not require any restriction on the growth of the cluster.

In this paper we show that there are a class of Rosenbluth like algorithms which yield a degeneracy of $`N`$ and which are straightforward to implement. The method provides an estimate of the number of lattice animals and can also yield estimates of any other desired properties of the animals such as their radius of gyration or perimeter multiplicities . We describe and justify the algorithm in Section 2 and present results to illustrate the use of the method in Section 3. Conclusions are given in Section 4

## 2 Algorithm

Any algorithm, suitable for the purpose of the enumeration of lattice animals using the Rosenbluth Monte Carlo approach, must satisfy two important criteria. First of all it has to be ergodic. That is to say, the algorithm should have a non zero probability of sampling any given cluster shape. The second criteria relates to the degeneracy that is associated with each cluster and requires this to be determinable. This degeneracy arises from the number of different ways that the same cluster shape can be constructed by the algorithm. While it is easy to devise methods of growing clusters that meet the first requirement, the second condition is more difficult to satisfy. For many simple algorithms the calculation of the degeneracy, for every cluster, can be a more complex problem than the original task of enumerating the number of lattice animals.

In the original Rosenbluth Monte Carlo approach of Care , this difficulty was overcome by ensuring that the degeneracy for all clusters of size $`N`$ was the same and equal to $`N!`$. However, to achieve this result the algorithm had to employ a somewhat elaborate procedure. This made the implementation of the method rather complicated, as well as limiting its possible extension to enumeration of other type of clusters. Here we shall consider an alternative algorithm, which while satisfying both of the above criteria, is considerably simpler than the algorithm proposed by Care. In Section 2.1 we describe the algorithm in its most basic form, before proving in Section 2.2 that the ergodicity and the degeneracy requirements are both met. In Section 2.3 we demonstrate how the basic algorithm can be further refined to improve its efficiency.

### 2.1 Basic Algorithm

Having chosen a suitable lattice on which the clusters are to be grown (square and simple cubic lattices were used in this study for 2D and 3D systems, respectively), a probability $`p`$ of acceptance and $`q=(1p)`$ of rejecting sites is specified. Although in principle any value of p between 0 and 1 can be selected, the efficiency of the sampling process is largely dependent on a careful choice of this value, as will be discussed later. In addition, an ordered list of all neighbours of a site on the lattice is made. For example, for a 2D square lattice this might read (right, down, left, up). While the order initially chosen is arbitrary, it is essential that this remains the same throughout a given run. In the basic algorithm, once chosen, the probability $`p`$ remains fixed during the Monte Carlo sampling procedure. However in Section 2.3 the effect of relaxing this requirement is discussed.

We construct an ensemble of $`N_E`$ clusters and for each of these calculate a weight factor which we subsequently use to calculate weighted averages of various cluster properties. For a property $`O`$ of the clusters, the weighted average is defined as

$$<O>_W=\frac{1}{N_E}\underset{\alpha =1}{\overset{N_E}{}}W_\alpha O_\alpha $$
(1)   
The weight associated with cluster $`\alpha `$ with $`N`$ sites is defined to be $`W_\alpha =1/(d_NP_\alpha )`$ where $`P_\alpha `$ is the normalised probability of growing the cluster and $`d_N`$ is a degeneracy equal to the number of ways of growing a particular cluster shape. It can be shown  that the weighted average can be used to estimate the number, $`c_N`$, of lattice animals of size $`N`$ and other properties such as the average radius of gyration $`\overline{R_N^2}`$:-

$`E[<1>_W]`$ $`=`$ $`c_N`$ (2)   
$`E[<R_\nu ^2>_W]`$ $`=`$ $`{\displaystyle \underset{\{\nu =1\}}{\overset{c_N}{}}}R_{N\nu }^2=c_N\overline{R_N^2}`$ (3)   
During the growth of each cluster we maintain a record of the sites which have been occupied, the sites which have been rejected and a ‘last-in-first-out stack’ of sites which is maintained according to the rules described below. Each cluster is grown as follows

1. Starting from an initial position, the neighbours of this site are examined one at a time according to the list specified above. An adjacent site is accepted with a probability p or else is rejected.
2. If the adjacent site is rejected, a note of this is made and the next neighbour in the list is considered.
3. If on the other hand it is accepted, then this becomes the current site and its position is added to top of a stack, as well as to a list of accepted sites. The examination of the sites is now resumed for the neighbours of this newly accepted site. Once again this is done in the strict order which was agreed at the start of the algorithm.
4. Sites that have already been accepted or rejected are no longer available for examination. Thus, if such a site is encountered, it is ignored and the examination is moved on to the next eligible neighbour in the list.
5. If at any stage the current site has no more neighbours left, that is all its adjacent sites are already accepted or rejected, then the current position is moved back by one to the previous location. This will be the position below the current one in the stack. The current position is removed from the top of the stack, though not from the list of accepted sites.
6. The algorithm stops for one of the following two reasons. If ever the number of accepted sites reaches N, then the algorithm is immediately terminated. In this case a cluster of size N is successfully produced. Note that unlike some of the other common cluster growth algorithms , it is not necessary here for every neighbour of the generated cluster to be rejected. Some of these might still be unexamined before the algorithm terminates. The second way in which the algorithm stops is when it fails to produce a cluster of size $`N`$. In this case, the number of accepted sites will be $`M<N`$, with all the neighbours of these $`M`$ sites already having been rejected, leaving no eligible sites left for further examination. From step (v), it is clear that in cases such as this, the current position would have returned to the starting location.
7. The probability of producing a cluster of size N, in a manner involving $`r`$ rejections, is simply $`p^{(N1)}q^r`$. Hence the weight, $`W_\alpha `$, associated with the growth of the cluster is given by

$$W_\alpha =1/(d_Np^{(1N)}(1p)^r)$$
(4)   
where the degeneracy, $`d_N`$, is shown below to be exactly $`N`$. Failed attempts have a zero weight associated with them. However they must be included in the weighted average of equation (1).
8. During the growth of a cluster of size $`N`$, we may also collect data for all the clusters of size $`M`$ where $`MN`$. It must be remembered that the weights for these smaller clusters must be calculated with a degeneracy of $`M`$.

A specific example is helpful in demonstrating the algorithm. Figure 1 displays a successful attempt in forming a cluster of size $`N=4`$, on a square lattice. The order in which the neighbours were examined was chosen to be right, down, left and up. Let us now consider various steps involved in construction of this cluster in detail. Beginning from the initial position labelled cell one, the adjacent site to the right of this position is examined. In this case the site is rejected and the current position remains on the cell one. Such rejected cells are indicated by the letter X. The next neighbour in the list is the one below, labelled cell two. As it happens this is accepted. Thus, the current position moves to this site and its position is added to the top of the stack, ahead of the position of cell one. The process of examining the neighbours is resumed for sites adjacent to cell two. Once again, following the strict order in the list, the site labelled three to the right of current position is considered first. This is also accepted and as before is placed at the top of the stack. At this stage the stack contains the positions of cells three, two and one, in that order. The current position is now cell three.

The site to the right of this, followed by the one below, are tested and both rejected in succession. Since both the neighbours to the left ( ie cell one) and the one above have already been considered, the current position has no more eligible neighbours left to test. Therefore, following step (v) above, site three is removed from the stack. This leaves the position of cell two at the top of the stack, making this the current position again. The cell two has two neighbours, the adjacent sites below and to the left, which are still unexamined. Of these, according to our agreed list, the site below takes precedent, but as shown in Figure 1 this is rejected. Current position remains on the cell two and the neighbouring site (cell labelled four) to the left of this position is tested. As it happens this is accepted. A cluster of the desired size $`N=4`$ is achieved, bringing this particular attempt to a successful end.

For the subsequent discussion, it is useful to represent a sequence of acceptance and rejections by a series of 1 and 0. Thus, for the case shown in Figure 1 we have {0,1,1,0,0,0,1}. Note that at any stage throughout a series, the position of the current site and that of the neighbour to be examined, relative to the starting cell, are entirely specified by the decisions that have been made so far. In other words, given a sequence of one and zeros we can determine precisely the shape of the cluster that was constructed. This is only possible because of the manner in which the neighbours of the current position are always tested in a strict pre-defined order. For an algorithm that considers the neighbouring sites at random, the same will clearly not be true.

The procedure described above needs to be repeated a large number of times, to obtain the weights for the ensemble average defined in equation (1). In particular, using equation 2, the number of lattice animals of size $`N`$ can now be determined.

### 2.2 Ergodicity and degeneracy of the algorithm

Let us now discuss the issue of the ergodicity of the algorithm. We wish to see whether, starting from any particular site on a given cluster, a series of acceptance and rejections (1 and 0) can always be determined which leads to that cluster shape. We stress that we are not concerned about how probable such a sequence is likely to be, but merely that it exists. We can attempt to construct such a sequence by following the same rules as our algorithm described above, with one exception; we accept and reject each examined site according to whether it forms part of the target cluster shape or not. Obviously, in the original algorithm, each such move has a non zero chance of occurring, provided $`p`$ is not set to zero or one. Since we only accept sites that belong to the cluster in question, it follows that if the sequence is successful then we would achieve the desired cluster shape. However, we might argue that for some choice of target cluster and starting position, a series started in this manner will always terminate prematurely. That is to say, it will inevitably lead to a failure, with only part of the required cluster having been constructed. Now, it is easy to see that this cannot be true. If the series fails, it implies that all the neighbouring sites of the sub-cluster formed so far are rejected. However, the rest of the cluster must be connected to this sub-cluster at some point. Hence, at very least, one neighbouring site of the sub-cluster must be part of the full cluster and could not have been rejected. Starting from any of the sites belonging to a cluster then, it is always possible to write down a sequence of one and zeros that will result in the formation of that cluster. Similarly, considering every starting point on a cluster of size $`N`$, another implication of the above result is that the corresponding cluster shape can be generated in a minimum of at least $`N`$ distinct ways.

Next, we shall show that the degeneracy of a cluster of size $`N`$ in our algorithm is in fact exactly $`N`$ (unlike the original algorithm of Care  which has a degeneracy of $`N!`$). Let us suppose that starting from a particular site on a given target cluster shape, our algorithm has two distinct ways of forming this cluster. Associated with each of these, a series of one and zeros can be written down, in the same manner as that indicated above. The two ways of constructing the cluster must necessarily begin to differ from each other at some stage along the sequence, where we will have a 1 in one case and a 0 in the other. Now since up to this point the two series are identical, the site being examined at this stage will be the same for both cases. This is rejected in one sequence (hence 0) whereas it is accepted in the other (hence 1). It immediately follows that these two differing ways of constructing the cluster cannot result in the same shape. Using this result, together with previous one regarding the ergodicity of the algorithm, we are lead to conclude that, starting from a given site on a cluster, the algorithm has one and only one way of constructing the cluster. Hence, for a cluster of size $`N`$, the degeneracy is simply $`N`$.

### 2.3 Refined algorithms

#### 2.3.1 Adjacent site stack

During the growth of the cluster a stack can be constructed of all the sites which are adjacent to the cluster and still available for growth. When a new site is added to the cluster, its neighbours are inspected in the predetermined sequence and any available ones are added to the top of this stack. (Note that this stack differs from that discussed in Section (2.1)). The choice of site to be occupied can be made from all the adjacent sites in a single Monte Carlo decision. Thus, if we consider the underlying process in the method described above, at each step there is a probability $`p`$ of the site being accepted and a probability $`q=1p`$ of the site being rejected. We therefore need to generate a random number with the same distribution as the number of attempts needed to obtain an acceptance. The probability of making $`k`$ attempts of which only the last is successful, is

$$p_k=q^{k1}p$$
(5)   
where $`1k<\mathrm{}`$ and $`_{k=1}^{\mathrm{}}p_k=1`$. In order to sample from this distribution we note that the associated cumulative distribution, $`C_m`$, is given by

$$C_m=\underset{k=1}{\overset{m}{}}q^{k1}(1q)=1q^m$$
(6)   
Hence if we generate a random number, $`\eta `$, uniformly distributed in the range $`0<\eta <1`$, then a number $`m`$ given by

$$m=\mathrm{Int}\left[\frac{\mathrm{ln}(\eta )}{\mathrm{ln}(q)}+1\right]$$
(7)   
will have been drawn from the required distribution. Thus we generate the number $`m`$ according to equation (7) and use this to determine which site on the stack is selected, with $`m=1`$ corresponding to the site at the top of the stack. If $`m>N_{adj}`$, where $`N_{adj}`$ is the number of available adjacent sites, the cluster growth is terminated as explained in step ((vi)) in Section 2.1. All the adjacent sites lying above the chosen site in the stack are transferred into the list of rejected sites. The list of adjacent sites is then adjusted to include the new available sites adjacent to the recently accepted site. As before, it is crucial that these are added to the top of the list in the strict predefined order.

#### 2.3.2 Variable probability

An apparent disadvantage of the methods so far described is that with fixed choice of probability, $`p`$, occasions arise when a cluster growth will terminate before reaching a cluster of size $`N`$, simply because the Monte Carlo choice rejected all the neighbouring sites. This problem can be overcome if the value of $`p`$ is allowed to vary as the cluster grows. The simplest method is to determine the number, $`N_{adj}`$, of available adjacent sites at each point in the cluster growth and select one of these sites with uniform probability. This effectively makes $`p=1/N_{adj}`$ and thereby increases the chances of growing a cluster of size $`N`$. Note that it is still possible for a cluster growth to become blocked. This happens when the chosen site is the one at the bottom of the current eligible neighbours list, thus causing all the other neighbouring sites in the list to be rejected in one step. If the newly accepted site has itself no unexamined neighbours to add to the list, the algorithm terminated prematurely. Modified in the manner described above the weight associated with a cluster is now

$$W_\alpha =\frac{\mathrm{\Pi }_{i=1}^NN_{adj}^i}{N}$$
(8)   
rather than the expression given in equation (4).

However, when this variable probability method was tested it was found that although it reduced the number of rejected clusters, it was inefficient at sampling the space of possible clusters when compared with method described in section (2.3.1). This inefficiency was measured by comparison of the standard deviation in the estimated cluster number for any given number of clusters in the sampling ensemble. It is thought that the inefficiency of the variable probability method arises because it gives too much weight to sites lower in the stack, yielding many non-representative clusters. It is possible that this problem could be overcome by using a non-uniform sampling distribution (cf ) but this was not tested in this work and the method described in (2.3.1) was used to obtain the results described in Section (3) .

## 3 Results

In order to test the algorithm described in Section (2) it was used to estimate the number of lattice animals on a square 2D lattice and a simple cubic 3D lattice for which exact results are known up to certain sizes . Before collecting data it was necessary to determine the optimum value of the probability $`p`$ with which an adjacent site is accepted during the cluster growth. The effect of changing $`p`$ on the estimated error in the number of clusters of size 50 on the 2D and 3D lattices can be seen in Figure 2. It can be seen that there is a fairly broad range of values of $`p`$ for which the error is a minimum and a value of $`p=0.6`$ was used to obtain the results described below for the 2D lattice and $`0.72`$ for the 3D lattice. The distribution of weights is log normal  and becomes highly skewed for large cluster sizes; this is a standard problem with Rosenbluth methods . The minimum in the error achieved by the choice of the value of the probability $`p`$ has the effect of minimising the variance of the distribution of the weights, $`W_\alpha `$.

In Table 1 we present results obtained using the algorithm defined in section 2 using the adjacent site stack method of section 2.3 to enumerate clusters on a simple cubic 3D lattice for clusters up to size 50. The results were obtained from an ensemble of $`2.5\times 10^7`$ clusters. The data took 3.3 hours to collect on a R5000 Silicon Graphics workstation using code written in the language C but with no attempt to optimise the code. Only $`30\%`$ of the clusters achieved a size of 50. The results are quoted together with a standard error, $`e^{est}`$, calculated by breaking the data into 50 blocks and determining the variance of the block means for each cluster size. If the number of samples in each block is sufficient, it follows from the central limit theorem that the sampling distribution of the means should become reasonably symmetrical. We therefore also quote a skewness, $`\xi `$, defined by 

$$\xi =m_3/m_2^{3/2}$$
(9)   
where $`m_i`$ is the $`i^{th}`$ moment about the mean of the sampling distribution. It is expected that $`\xi \stackrel{<}{}0.5`$ for a symmetrical distribution and $`\xi >1`$ for a highly skew distribution. The statistic $`\xi `$ should be treated with some caution since it is likely to be subject to considerable error because it involves the calculation of a third moment from a limited number of data points.

Exact results are known for clusters up to size 13  and in the table we quote the values for the quantity $`\chi `$ defined by

$$\chi _M=|\frac{c_M^{exact}c_M^{est}}{c_M^{exact}e_M^{est}}|$$
(10)   
where $`c_M`$ is the number of clusters of size $`M`$ and it can be seen that all the values of $`\chi `$ are $`O(1)`$. Hence we assume that $`e^{est}`$ is an acceptable method of estimating the error in the method. However it is likely that the $`e^{est}`$ will underestimate the true error if the distribution becomes more skew. We also quote in Table 1 the values of $`c_N`$ calculated by Lam  using a Monte Carlo incomplete enumeration method together with the error estimates reported for this method.

In Table 2 we quote data collected from a square two dimensional lattice by collecting data from $`2.5\times 10^7`$ clusters up to size 50. This data only took 1.45 hours to collect but only $`2\%`$ of the clusters achieved a size of 50. Comparison is given with exact results  up to clusters of size 19. The rate of growth of errors for the two and three dimensional data is shown in Figure 3 and it can be seen that the errors associated with the method diverge are beginning to diverge quite rapidly above clusters of size 50. This behaviour is to be expected with a technique which is based on sampling from a log normal distribution. In the previous paper  equivalent results were obtained for clusters up to size 30 with approximately the same sample size. The improvement up to clusters of size 50 obtained by the new method arises because the weight associated with clusters of a certain size is generated from roughly half as many random numbers. This effectively halves the standard deviation of the log normal distribution of the weights and allows larger clusters to be sampled before the method becomes unusable.

## 4 Conclusions

We have described a simple Rosenbluth algorithm for the Monte Carlo enumeration of lattice animals and clusters which can be applied to any lattice topology. A merit of the scheme is that for thermal systems it may be easily adapted to include Boltzmann weightings following, for example, the arguments used by Siepmann at al  in the development of the configurational bias technique. Similarly, the method can be applied to calculation of the averaged properties of a cluster of a given size, in the site percolation problem. In this case we have

$$<O>=\frac{<(1P)^tO>_W}{<(1P)^t>_W}=\frac{_{\alpha =1}^{N_E}W_\alpha (1P)^{t_\alpha }O_\alpha }{_{\alpha =1}^{N_E}W_\alpha (1P)^{t_\alpha }}$$
(11)   
where $`P`$ is the probability of site occupation in the percolation problem of interest and $`t_\alpha `$ the number of perimeter sites  of the cluster $`\alpha `$. Preliminary results also indicate that the method may be useful in the study of the adsorption of clusters onto solid surfaces. A possible numerical limitation of the method arises from the highly skew probability distribution of Rosenbluth weights which occurs for large cluster sizes. However the method presented in this work is able to work to considerably higher cluster sizes than the one described in  before this becomes a problem.