Application of the parallel multicanonical method to lattice gas condensation

We present the speedup from a novel parallel implementation of the multicanonical method on the example of a lattice gas in two and three dimensions. In this approach, all cores perform independent equilibrium runs with identical weights, collecting their sampled histograms after each iteration in order to estimate consecutive weights. The weights are then redistributed to all cores. These steps are repeated until the weights are converged. This procedure benefits from a minimum of communication while distributing the necessary amount of statistics efficiently. Using this method allows us to study a broad temperature range for a variety of large and complex systems. Here, a gas is modeled as particles on the lattice, which interact only with their nearest neighbors. For a fixed density this model is equivalent to the Ising model with fixed magnetization. We compare our results to an analytic prediction for equilibrium droplet formation, confirming that a single macroscopic droplet forms only above a critical density.


Introduction
Despite the large amount of existing literature, condensation remains under investigation up to today by means of analytical treatment, numerical studies, and experiments [1,2,3,4,5,6,7,8]. In addition to the increasing interest in non-equilibrium properties and droplet size distributions, questions regarding the equilibrium properties remain open as well. For example, Biskup et al. [3,4] showed that in equilibrium an over-saturated gas system is either in the evaporated phase or in the condensed phase consisting of a single macroscopic droplet. The probability for droplets of intermediate size is negligibly small. This result was proven for the two-dimensional Ising model and numerically verified [7,8]. We set out to test this analytical result also for the three-dimensional case of a lattice gas. Due to the nucleation barrier, we applied a parallel implementation [9] (for related methods see [10,11,12]) of multicanonical simulations [13,14,15] to circumvent a possible hysteresis, driving the system back and forth between the condensed and evaporated phase.
We start by describing the main method and the theoretical arguments together with the model. This is followed by a comparison to micromagnetic simulations of the equivalent Ising model in order to verify the correctness. We discuss the speedup with increasing parallelization and compare the simulation effort with the independent micromagnetic simulation. Afterwards, the results in two and three dimensions are presented and discussed, followed by some concluding remarks.

Method
The multicanonical method [13,14,15] enables the sampling of a broad parameter space in a single simulation by modifying the acceptance probability according to a simple idea. The Boltzmann factor is replaced, or sometimes extended, by a weight function with the purpose to increase the probability to reach states that are suppressed otherwise. This can be achieved for a variety of order parameters leading to multimagnetic, multibondic and many more realizations but can best be explained for the case of the internal energy E. Formally, this can be described by writing the partition function in terms of the density of states Ω(E), where we have already transformed the Boltzmann factor into the weight function W (E). Analogous to the Metropolis algorithm, new states (E n ) are generated from old states (E o ) and accepted with If we were to know Ω(E), we could sample the full energy range with a flat distribution choosing W (E) = 1/Ω(E). However, this is in general not the case such that W (E) has to be obtained in an iterative way. This can be simple or complex; it can modify the weights on the fly or rely on equilibrium distributions in each iteration. A simple example would be to perform equilibrium simulations with initial weights W i and estimate the consecutive weights . Hence, if a state occurs more often the weight is reduced more strongly.
We applied a parallel version of this multicanonical method which has been shown to exhibit an ideal scaling for the Ising model and a very good scaling for the q = 8 Potts model [9]. This parallel implementation relies on equilibrium iterations, distributing the sampling of the distribution H i (E) on p independent processes. Since the weights are not modified within a single iteration, each process returns a contribution to the same distribution and they can simply be summed up: Afterwards, the consecutive weights are calculated from the previous weights together with the total histogram according to whatever rule one may choose, in our case by a recursive scheme [15]. Altogether, this leads to an efficient parallelization due to a very limited amount of communication. Below, we will show the applicability of the method to a general problem, here the condensation of a lattice gas in two and three dimensions. As in real life and other than in [9], we will not pay attention to the optimization of each degree of parallelization but want to demonstrate the value of this parallelization to modern computational statistical physics.

Condensation of a lattice gas
A general theory of equilibrium droplet condensation was proposed by Biskup et al. [3,4], with a rigorous analytic solution for the two-dimensional Ising spin model. In principle they describe the interplay of entropy maximization by vacuum fluctuations with energy minimization by forming a droplet. They showed that the probability for intermediate-size droplets vanishes with increasing system size. The free energy consists of a contribution from the vacuum fluctuation of excess particles δN with the isothermal compressibilityκ = βκ = β (N − N ) 2 /V and a contribution from the single largest droplet of size where τ W is the surface free energy of a (Wulff shaped) droplet of unit volume (see Fig. 1). This was translated into the language of the Ising spin model in [3,7]. Here, we will concentrate on the formulation of a lattice gas, where each lattice site is either occupied by a particle or empty, i.e. n i ∈ {0, 1}. For a fixed number of particles this model is equivalent to the Ising model with spins s i ∈ {−1, 1}, at fixed total magnetization, when s i = 2n i − 1: Thus, for fixed coupling constants, a scaling of the temperature by T Ising = T I = 4T and the addition of a system-size dependent constant transforms the energy scale of an Ising model into the corresponding one of a lattice gas. All other observables are then mapped accordingly by scaling the temperature. The contributions to the free energy can as well be formulated in the language of a lattice gas. Consider a temperature-dependent vacuum or background density ρ 0 = N 0 /V , which may be identified with the spontaneous magnetization in the Ising model via Adding particles to the system leads to a total particle excess δN = N − N 0 . According to Biskup et al. this excess δN may be decomposed into the particle excess inside the droplet δN D and the remaining particle excess in the fluctuating phase δN F . For the infinite-size system and from the equivalence to the Ising model follows that the gas density is given by ρ G = ρ 0 and the liquid density by ρ L = 1 − ρ 0 . Hence, in a droplet of size V D we expect an excess of particles δN D = (ρ L − ρ G )V D . Introducing the fraction of excess particles inside a droplet λ = δN D /δN , allows us to rewrite the decomposed particle excess as fractions δN D = λδN and δN F = (1 − λ)δN . Altogether the total free energy F = F drop + F fluc then reads Rewriting the particle excess in terms of a volume and with the dimensionless parameter Whileκ has the same value as χ in the Ising model, τ W has to be rescaled such that τ I W = 4τ W . At fixed temperature, the parameters ρ 0 ,κ, τ W are constant. Minimizing eq. (9) leads to the infinite-size solution for the equilibrium volume fraction λ = λ(∆) (see [3,4]). In order to compare to the analytic predictions, we have to measure the average size of the largest droplet at different densities but fixed temperature. Because of the increasing correlation times in large systems, we applied the multicanonical method, driving the system between the evaporated and condensed phase, and reweighted to the relevant temperatures afterwards. The energy ranges were obtained through short parallel tempering [16] runs in the beginning of each simulation.

Results
In order to assess the parallel algorithm, we performed two independent simulations of equivalent but different models, namely the Ising spin model and the lattice gas. While the lattice gas was simulated with the parallel multicanonical algorithm explained above, the Ising model was simulated with a Metropolis algorithm at fixed magnetization, also referred to as micromagnetic simulation. Both approaches needed independent simulation runs for each measurement point. We chose the inverse temperature in the Ising language to be β I = 0.369, which corresponds to a temperature T = 0.6995 in the particle picture. Figure 2 shows very good agreement between the two completely independent and different methods.
A direct comparison of simulation times is rather difficult. For one, the multicanonical method provides additional information for various temperatures, and this comes at a price: A time-consuming weight iteration, which is not contributing to the production run. Whether or not the overhead of the initial weight iteration is justified, will depend on the desired size of the statistical error. Furthermore, the approaches will scale differently with system size, which we  Comparison of simulation data from micromagnetic Metropolis simulations of the Ising model and parallel multicanonical simulations of the lattice gas in three dimension. Here, the volume is fixed at L = 30 while the number of spins/particles are modified.
do not want to address at this point. For the plot in Fig. 2, the average simulation time per data point of both methods was about half an hour on common 2.5GHz cores, but the multicanonical simulation used 36 cores instead of one in the case of the micromagnetic simulation. However, with increasing system size the free-energy barrier between the evaporated and the condensed phase increases [17], and the necessary computation time for the single-core micromagnetic simulation exceeds the limits of practicability.
Next, we discuss the speedup of the parallel method with increasing degree of parallelization. As mentioned above, we intentionally did not pay attention to optimization of the simulation parameters. Rather we want to demonstrate the applicability of the parallelization by choosing appropriate parameters from single-process simulations. Hence, we fixed the total number of sweeps per iteration to 25600 such that each process performed 25600/p sweeps. A single sweep consists of L d updates. The number of particles was chosen to be N = 1000 and the energy range adapted accordingly. The simulation was defined to be converged if a total of 30 "tunneling events" were counted, meaning that the simulations had traveled at least 30 times from the maximal energy to the minimal energy or vice versa.
When we want to estimate the speedup, we need to consider statistical averages. This is due to the fact that the weight iteration relies on the sampled data such that the number of iterations until convergence cannot be predicted and may, moreover, vary depending on the initial seed of the random number generator. To this end, we performed 32 simulations for each degree of parallelization. The effect is best seen in the average number of iterations until convergence ( Fig. 3 (a)). An increasing number of independent Markov chains improves each estimate of the distribution belonging to the current weights. This optimizes the estimation for the successive weights and reduces the total amount of necessary number of iterations. Of course, this is in general not true on all scales [9], for example when the number of sweeps per iteration becomes too small or if the optimal number of sweeps is more complex than the assumed 1/p behavior.
The speedup is defined in terms of the average time until convergence; it is given by the ratio of a single-core to a p-core simulation: For a fair comparison, both times were obtained with the same program. Because the communication between processes occurs only at the end of every iteration, the speedup should be ideal with a linear behavior of slope one, i.e., doubling the amount of processes should speed up the time by a factor 2. We can see in Fig. 3 (b) that this is indeed true for the case of three dimensions and even better for two dimensions. This may be explained by the reduced number of iterations, due to the independence of Markov chains.
Having demonstrated the correctness of our model and the efficiency of our method, we want to apply it to the problem of condensation in a lattice gas. Figure 4 shows the results for selected temperatures in two and tree dimensions with up to 24 cores in each simulation. The temperature in two dimension was chosen to correspond to the equivalent Ising temperature T I = 1.5 used in Ref. [7], allowing us to use the same parameters. In three dimension we chose a temperature T I = 2 below the roughening transition T I R ≈ 2.4537 [18], where lowtemperature approximations hold. This allowed us to use a low-temperature series expansion for the interface tension: σ = 1.9072 . . . [19]. Furthermore, below the roughening transition the condensate tends to form a cube, such that the surface free energy may be assumed as τ I W ≈ 6σ. This is the result for the Ising model and has to be rescaled by the same factor as the temperature, hence τ W ≈ 6σ/4 (the same argument applies to the two-dimensional case). The remaining free parameters, ρ 0 = (1 − m 0 )/2 andκ = χ, were obtained by standard Metropolis simulations of the Ising model on very large lattices. An overview of the involved parameters is given in Table 1.
Opposite to Ref. [7], we always fixed the number of particles and varied the density by modifying the volume of the system. We can see in Fig. 4 that this choice leads to a similar behavior in two dimensions as in Ref. [7], where the predicted ∆ c separates the evaporated from the condensed phase quite well. Of course, we see finite-size effects as reported before, but the qualitative picture is satisfying. In three dimensions, on the other hand, we observe a strong deviation of small systems from the predicted critical value of the parameter ∆. This becomes weaker but is still present for the largest systems investigated. Specifically, the system needs to go to larger ∆ or, equivalently, particle density in order to condensate. With increasing system Table 1. Converted parameters for the lattice gas system in two [8] and three dimensions, where for the latter case ρ 0 = (1 − m 0 )/2 andκ = χ are obtained from Metropolis simulations of the Ising model at various system sizes and τ W ≈ 6σ/4 with σ from a low-temperature expansion [19].  Condensation of a lattice gas in two (left) and three (right) dimensions. The temperature is chosen, such that the two-dimensional case corresponds to the Ising model at T I = 1.5 in Refs. [7,8] and the three-dimensional case is below the roughening transition. The differently colored curves belong to fixed number of particles N , while the density is varied by changing L. Each data point corresponds to an independent simulation with up to 24 cores. size, the critical value of ∆ decreases, while the shape of the curves seems to converge to the predicted functional shape. This may be a finite-size effect that is more pronounced in three dimensions, which calls for more detailed studies in future work. In both cases, we do not have any free parameters left.

Conclusion
We have demonstrated the applicability of an efficient but simple parallelization of the multicanonical method to the case of lattice gas condensation in two and three dimensions. In these cases, the speedup was shown to scale ideally (or even better) due to the fact that independent Markov chains improve the sampling of equilibrium distributions.
Applying this method to the problem of condensation, we could reproduce the results from [7] in two dimensions and present new results in three dimensions which indicate large finitesize effects. In the latter case, condensation systematically occurred at densities higher than predicted. While the deviation decreases with system size, the available data does not allow for quantitative predictions due to the low resolution at the critical density caused by varying the volume. This has still to be investigated in more detail.