A continuous fractional component (CFC) approach increases the probability of particle swaps in the context of vapor-liquid equilibrium simulations using the Gibbs ensemble Monte Carlo algorithm. Two variants of the CFC approach are compared for simulations of pure Lennard-Jones (LJ) fluids and binary LJ mixtures as examples. The details of an exemplary CFC implementation are presented. Recommendations are provided to reduce the effort required for the suggested problems.

## I. INTRODUCTION

Since its original development in the 1980s, the Gibbs ensemble Monte Carlo (GEMC) algorithm has been successfully applied to the simulation of fluid phase equilibria. Among the most important advantages of the GEMC algorithm are its simple, intuitive methodology, and the fact that data for both coexisting phases can be obtained from a single simulation. The data include macroscopic properties such as the phase molar fractions, *x _{i}*, the phase densities,

*ρ*, and the chemical potentials,

_{i}*μ*.

_{i}^{1}A broad range of open-source GEMC programs has been developed.

^{2–6}However, a significant limitation of the GEMC algorithm is that the acceptance probability for molecule swaps becomes low for insertions into dense phases, especially if components with large differences of molecular sizes are simulated.

The goal of this paper is to discuss several variants of the GEMC algorithm, which increase the probability of particle swaps in the context of vapor-liquid equilibrium simulations. For this purpose, two variants of the continuous fractional component (CFC) approach are compared for pure Lennard-Jones fluids, binary Lennard-Jones mixtures, ethane, and carbon dioxide. Details of a continuous fractional component program in mathematica are presented. Simulation protocols are recommended to reduce the effort required for the proposed problems.

This paper is structured as follows: In Sec. II, the principles of the Gibbs ensemble Monte Carlo approach are presented. The challenge of inserting molecules into dense phases is discussed in Sec. III, and two frequently used variants of the continuous fractional component approach are introduced in Sec. IV. The implementation of these algorithms is presented in Sec. V. In Sec. VI, the continuous fractional component algorithm is compared to the conventional GEMC approach by discussing the acceptance probabilities, followed by a comparison of GEMC and the continuous fractional component algorithm and an assessment of the implementations regarding stability, comparison of simulation results to the literature, and the computational performance. Conclusions are given in Sec. VII, and Sec. VIII provides suggested problems for interested readers.

## II. GIBBS ENSEMBLE MONTE CARLO

*r*is the position of particle

_{j}*j*, and the summation is over all molecules

*N*in the system. For spherically symmetric molecules, the Lennard-Jones (LJ) potential function (cf. Fig. 3) can be used.

^{1}The LJ-parameters can be determined by ab initio algorithms

^{8}or from experimental data.

^{9}

*ξ*, which is uniformly distributed on (0, 1). If

_{i}*ξ*is greater than the calculated acceptance probability, the move is rejected.

_{i}In some implementations of the acceptance algorithm, an early rejection method called “overlap” is used, where an overlap radius $ r ov$ is defined. If the distance between two interaction sites is smaller than $ r ov$, the trial move is automatically rejected before applying Eq. (4).^{10}

A certain number of trial moves is executed during the three phases of the simulation: (i) in the warm-up phase, only displacements are attempted to “melt” the initial arrangement of the molecules, which usually are on a grid; (ii) in the equilibration phase, the system is equilibrated using the three trial moves; and (iii) during the production phase, a sufficient number of states is sampled to obtain ensemble averages of the macroscopic properties of the boxes such as the phase molar fractions *x _{i}*, the phase densities

*ρ*, and the chemical potentials

_{i}*μ*. Samples are taken only after a certain amount of trial moves have been executed. The routine of executing a set number of trial moves and sampling the system once is referred to as a Monte Carlo cycle.

_{i}^{10}

## III. THE CHALLENGE OF INSERTING MOLECULES INTO DENSE PHASES

Practical applications of the usual GEMC algorithm are limited by the low acceptance probability of molecule swaps for insertions from the vapor box into the liquid box, especially if molecules with large differences of their molecular sizes are simulated. The lower the temperature and the higher the density, the greater the number of swaps that need to be attempted, because the probability of accepting a swap decreases. This low swap acceptance probability is mainly due to the strong repulsive force when two particles become very close (cf. Fig. 3), and a randomly chosen insertion position is increasingly likely to be close to another molecule with increasing density. As an appropriate number of attempted swaps per Monte Carlo (MC) cycle, it has been proposed to start with half the number of molecules in the system and determine an appropriate number in iterative pre-simulations. For example, if the system contains 500 molecules, 250 swap moves are proposed as the initial number of swap moves per cycle.^{11} If the chemical potential is insufficiently equilibrated because not enough swap moves were accepted by the end of the simulation, the number of attempted swap moves is increased before running the simulation again.

However, at low temperatures and accordingly low acceptance probabilities, simply increasing the number of attempted swap moves is insufficient, because each trial move comes at a computational cost causing the duration of a simulation to become excessively long.

Several approaches have been proposed to overcome this problem. One of the earliest improvements is based on the idea that the acceptance probability for swaps of small particles is significantly higher than that of larger ones. Therefore, only the smaller particles are swapped between the boxes. The larger ones are not swapped directly but exchanged with the smaller molecules. As displayed in Fig. 2, this exchange is done by changing the identity of a small particle into a large particle in one box and the other way around in the second box. The results of this approach are in good agreement with experimental data for systems of noble gases such as argon, krypton, neon, and xenon.^{12}

The “inflating-flea” algorithm is similar to this approach. In the inflating-flea algorithm, a few small particles, which can be easily swapped between the boxes, are added. These “entrainer particles” are then swapped with the actual particles in a separate step. The main limitation of this algorithm is that the system is no longer identical to the original system.^{13,14}

In the configurational bias GEMC algorithms, the random insertion of an entire molecule is replaced by the gradual insertion of its constituting atoms, one after another in a way that avoids unfavorable, e.g., overlapping, positions. This insertion is done in such a way that preferable configurations, which are more likely to be accepted, are sampled. Configurational bias GEMC algorithms have been successfully applied to alkane mixtures.^{15,16} A substantially different approach, which avoids particle swaps, combines the GEMC approach and thermodynamic integration of the Gibbs–Duhem equation along the saturation line.^{17}

## IV. THE CONTINUOUS FRACTIONAL COMPONENT ALGORITHMS

Methodologically, the continuous fractional component algorithm is similar to the inflating-flea algorithm. Instead of adding small tracer particles, continuous fractional component algorithms are based on the introduction of fractional molecules. The interactions of a fractional with the rest of the system are scaled, so that the particle can be inflated or contracted gradually. By doing so, a gradual insertion or deletion replaces the conventional GEMC swap move for which an entire molecule is moved from one box to the other at once. In this way, a reasonable acceptance probability of particle swaps into the dense phases is achieved. A comprehensive overview of the development of the continuous fractional component algorithm can be found in Ref. 18. We summarize the basic ideas of its methodology in the following.

*λ*is introduced that controls the interactions of the fractional particle. The interaction of the fractionals with the other particles is scaled as

*ξ*is an adjustable parameter, which is usually set to 0.5,

^{19}

*r*is the distance between the Lennard-Jones particles, and

_{ij}*σ*and

_{ij}*ϵ*are the size and energy pair parameters. They are calculated with an appropriate combining rule just like for a conventional GEMC simulation, for example, with the Lorentz-Berthelot rules given in Eqs. (6) and (7). The scaled Lennard-Jones potential is illustrated in Fig. 3.

_{ij}^{20}

As can be seen in Fig. 3, the interaction of a fractional with another particle is small when *λ* is small. Thus, the insertion probability of the fractionals into a dense phase is high, because the energy contribution, even for very close particles, is low.

Not all fractionals have the same value of *λ*. The *λ* distribution can become steep because states with low *λ* values are unlikely. We use the Wang–Landau approach to determine a biasing function $ \eta ( \lambda )$ so that all *λ* values are equally likely.^{21,22} This biasing function $ \eta ( \lambda )$ is initialized to zero at the start of the simulation and is determined iteratively during the equilibration phase. For this purpose, a biasing increment *v* is used, and the values of *λ* are divided into ten equal sized bins, e.g., the first bin is $ 0 < \lambda < 0.1$. Each time a new *λ* bin is visited, the respective bin and, thus, $ \eta ( \lambda )$ is incremented by $ ln \u2009 ( v )$. For example, if a trial move resulting in $ \lambda = 0.15$ is accepted, the value of the second bin is incremented by $ ln \u2009 ( v )$. To ensure a smooth *λ* distribution, a linear interpolation is used between the bins. In this way, trial moves resulting in *λ* bins that have been visited often are less likely to be accepted, and the *λ* distribution becomes flatter.^{23}

*v*is adjusted. If each bin, which should on average contain 10% of the visited states, contains at least 1% of all

*λ*states that have been visited during the respective evaluation period, the biasing increment

*v*is decreased according to

^{21,23,24}

### A. Continuous fractional component with two fractionals

First, the proposed change of the coupling value $ \Delta \lambda $ is randomly sampled in the range $ [ 0 , \Delta \lambda max ]$, where $ \Delta \lambda max$ is chosen using the algorithm by Allen and Tildesley^{1} to create a 50% continuous fractional component acceptance probability.

Separate biasing functions $ \eta ( \lambda )$ are used for each box to guarantee a flat *λ* distribution in each box. Once a change in *λ* has been sampled, we check to see if the new value of *λ* is within the range [0, 1]. If so, we attempt to change the coupling parameter *λ* of both fractional molecules in their respective boxes. Both moves are illustrated in Fig. 4.

*λ*is outside the range $ [ 0 , 1 ]$, a molecule swap is attempted as illustrated in Fig. 4(b). The new values for

*λ*of the fractionals are calculated using

*λ*close to unity is small. At the same time, the negative energy contribution from removing the fractional with the small coupling parameter is also small and, most importantly, the interactions of the newly inserted molecule are minimal. Consequently, the acceptance probability of this procedure of exchanging a particle to the other box is high and favorable to facilitate particle swaps. The change of the interaction energies

*ν*and the biasing values

*η*are used to calculate the acceptance probability given by

^{19}

### B. Continuous fractional component with one fractional

A continuous fractional component algorithm that is similar to CFC2F and is referred to as CFC1F uses only one fractional per component instead of two coupled ones.^{24} Consequently, the fractional is present in only one of the two boxes. A third continuous fractional component trial move—the “exchange move”—is introduced that swaps the fractional from one box to the other. Another difference is that the type of the continuous fractional component trial move is randomly chosen and independent of the current value of *λ*: 50% of the continuous fractional component trial moves are *λ* changes, 25% are exchanges, and 25% are swaps.

*λ*is similar to that of CFC2F. The main difference is that only the change of the interactions of one particle needs to be calculated. The acceptance rule is

^{24}

*λ*and the positions of the particles remain unchanged, but the number of full particles in the boxes can change, is

*i*in the old and new boxes. This trial move is likely to be accepted for values of

*λ*close to unity, where the fractionals' interactions are almost equal to those of a full particle.

^{24}The continuous fractional component exchange trial move is displayed in Fig. 5(a).

*λ*and the number of full molecules remain constant, but the fractional is moved from one box to the other, the acceptance rule is

*λ*for which the interactions of the fractional are weak.

^{24}

## V. IMPLEMENTATION OF THE ALGORITHM

The simulations were conducted using software optimized for comparing the two continuous fractional component algorithms we have discussed. The software is based on the GEMC algorithm proposed by Panagiotopoulos,^{10,25,26} comprising a warm-up, equilibration, and production phase.^{1} Analytical tail corrections account for the truncation of the intermolecular interactions of the Lennard-Jones potential.^{27} The pressure is determined from the ideal and viral contributions.^{28} The Widom ghost particle algorithm is used to determine the chemical potential.^{29} For the calculation of cross-interactions, the geometric mean or the Lorenz-Berthelot combining rules [see Eqs. (6) and (7)] are available.^{30} For electrostatic long range corrections, the Ewald or Wolf summations are implemented.^{31,32} This simulation package, written in mathematica,^{33} is part of the supplementary material^{34} and also available in Ref. 35.

## VI. ASSESSMENT

To have a reasonable computation time and to reduce finite-size effects, periodic boundary conditions are used. All the simulations were conducted with 500 molecules, using the Lennard-Jones parameters for argon: $ \sigma = 3.3952 \u2009 \xc5 , \u2009 \u03f5 / k B = 116.79 \u2009 K$.

Additionally, three binary Lennard-Jones mixtures ranging from $ \sigma 22 / \sigma 11 = [ 1 \u2212 0.5 ]$ and $ \u03f5 22 / \u03f5 11 = [ 0.75 \u2212 0.5 ]$ were studied at the reduced temperatures $ T * = k B T / \u03f5 = 1.0$ and 0.75, where *σ _{ii}* and

*ϵ*designate the Lennard-Jones size and energy parameters of the pure components.

_{ii}^{36}The Lennard-Jones parameters of the first component were set to unity ( $ \sigma 11 = 1 \u2009 \xc5 , \u2009 \u03f5 11 / k B = 1 \u2009 K$), and Lorentz-Berthelot combining rules were used. The Lennard-Jones fluids represent simple molecules such as noble gases satisfactorily, because such components cannot deform and, thus, have no intra-molecular energy.

More complex molecules were also studied, including the two-center molecule ethane, both as a Lennard-Jones-particle and as a rigid two-center molecule. Furthermore, carbon dioxide was simulated as an example for partially charged molecules. The parameters used for these molecules are provided in the supplementary material,^{34} and the molecular models are visualized in Fig. 6.

_{2}are provided in Ref. 34. For the quantitative assessment, reduced LJ units, which are often used in the literature, are used [see Eqs. (17)–(20)]. For mixtures, the size parameter

*σ*is calculated using Eq. (21)

### A. Continuous fractional component acceptance probabilities

To investigate the acceptance probabilities of the more complex continuous fractional component algorithms, we applied the CFC1F method for argon as a pure Lennard-Jones-fluid at $ T * = 0.7$. After each Metropolis check, we recorded if the trial move was accepted, along with the value of *λ* of the fractional and the distance to the closest particle. With these data, the cumulative distribution of finding an insertion position within a certain distance of the closest molecule was calculated for the vapor and liquid boxes, as illustrated in Fig. 7(a), showing the distributions in the liquid and the vapor phase on different horizontal axes.

As expected, the distance between the molecules in the vapor phase is significantly larger than that of the liquid phase. The distribution in the vapor phase is broad compared to that of the liquid phase and covers a distance of $ 1 \u2272 r / \sigma \u2272 8$. In the liquid phase, the distribution is more narrow and covers only a range of $ 0 \u2272 r \sigma \u2212 1 \u2272 0.8$, which is reasonable given that $ T *$ is low.

The acceptance probabilities for insertions are shown for the liquid phase in Fig. 7(b). Based on Fig. 7(a), the distance to the nearest neighbor was discretized into four equidistant bins in the range $ 0 < r / \sigma < 1$; for example, the first bin covered the distance $ 0 < r / \sigma < 0.25$. The values of *λ* were also discretized into four equidistant bins, the first bin covering $ 0 < \lambda < 0.25$. We observed that, for $ \lambda < 0.25$ and the nearest neighbor closer than $ 0.75 \u2009 \sigma $, the acceptance probability is $ \u2248 0.5 %$. In contrast, not one attempted insertion was accepted when $ \lambda > 0.25$ at the same distance to the nearest neighbor. Swaps with $ \lambda > 0.25$ were only accepted when randomly occurring cavities were hit, as indicated by the large distance to the closest neighbor ( $ 0.75 < r / \sigma < 1$). This low acceptance when two molecules are very close is due to the strong repulsive forces resulting from the $ ( \sigma / r ) 12$ part of the Lennard-Jones potential.

Attempted exchanges with $ 0.75 < \lambda < 1$, which is similar to the non-adjusted potential function (cf. Fig. 3), were also only accepted when the distance to the closest neighbor was greater than $ 0.75 \u2009 \sigma $. Thus, conventional GEMC swaps are successful only if very rarely occurring cavities are hit. Consequently, many insertion positions have to be sampled and tested to facilitate swaps.

In contrast, insertions even into very dense phases are accepted more often using the CFC1F algorithm, and the strong repulsive forces at close distances of the Lennard-Jones potential do not inhibit the acceptance probability of swaps as much.

### B. Computational performance

Given that the goal is to facilitate swaps more efficiently, the computation time needed to perform such moves and the acceptance probabilities of the algorithms are important. In Table I, the computational times needed for the different parts (overlap check, Metropolis check, update of the variables—see Sec. V) of the different trial moves are compared to the time needed to perform one displacement. This normalization is especially useful if simulations on different computers are compared because the computational time depends on the hardware and software of the machine used.

Method . | Type . | Overlap . | Metropolis . | Complete . |
---|---|---|---|---|

GEMC | Translation | 0.5 | 0.6 | 1.0 |

$ \Delta V$ | n.a. | 72.0 | 74.4 | |

Swap | 0.5 | 0.7 | 1.8 | |

CFC 2F | Lambda | n.a. | 1.1 | 2.4 |

Swap | 0.4 | 1.9 | 4.9 | |

CFC 1 F | Lambda | n.a. | 0.4 | 0.9 |

Exchange | n.a. | 1.5 | 4.7 | |

Swap | 0.5 | 0.6 | 1.0 |

Method . | Type . | Overlap . | Metropolis . | Complete . |
---|---|---|---|---|

GEMC | Translation | 0.5 | 0.6 | 1.0 |

$ \Delta V$ | n.a. | 72.0 | 74.4 | |

Swap | 0.5 | 0.7 | 1.8 | |

CFC 2F | Lambda | n.a. | 1.1 | 2.4 |

Swap | 0.4 | 1.9 | 4.9 | |

CFC 1 F | Lambda | n.a. | 0.4 | 0.9 |

Exchange | n.a. | 1.5 | 4.7 | |

Swap | 0.5 | 0.6 | 1.0 |

The conventional GEMC moves efficiently utilize the overlap optimization for displacements and swaps moves: For translations, approximately 20% of the computation time can be saved when an overlap is detected. Apart from that, a substantial time—approximately 40%—is needed to update the state vectors, where the information of the system (e.g., the positions of the individual atoms and molecuels are stored), and intermolecular energies when a proposed move is accepted. The reason is that the interactions of the moved molecule need to be updated in both the box where it is deleted and the box where it is inserted. Consequently, an accepted GEMC swap consumes almost twice as much computation time than a translation.

The volume change is the slowest GEMC move, because the interactions of all the molecules need to be calculated and updated.

The CFC2F algorithm is approximately two times slower compared to the conventional GEMC swap move due to the coupling of the fractionals, because the interactions of both fractionals need to be evaluated for every proposed *λ* change. This effect is especially pronounced when a swap move is accepted, because then the change of the interaction energies of four molecules must be calculated and updated.

In comparison, the CFC1F algorithm involves less computational effort. Moves adjusting the coupling strength of the fractional are faster compared to CFC2F moves because, just as for a displacement, only interactions of one molecule in one box need to be calculated. Exchange moves, however, are slower because the interaction of the new full molecule needs to be taken into account as well. Swap moves are more or less equivalent to the GEMC swap moves, and no significant difference in the computation time is observed.

### C. Argon as a LJ fluid

In the following, we compare the GEMC and the two continuous fractional component algorithms using three criteria: (i) The stability of the simulations as characterized by the standard deviation of the number of molecules per box; (ii) the agreement of the simulation results with published results; and (iii) the total time for one successful swap. In contrast to Table I, where only the computational time needed for performing one accepted swap is listed, Tables II and III show the total time, including the time needed for all non-accepted swaps and the time needed to update the state variables when the move is accepted.

. | . | . | $ t Swap$ . | $ \Delta \rho V *$ . | . | $ \Delta P V *$ . | . | $ std ( n V )$ . |
---|---|---|---|---|---|---|---|---|

System . | $ T *$ . | Algorithm . | 1 . | $ \sigma 3 10 \u2212 3$ . | $ \Delta \rho L *$ . | $ \sigma 3 \u03f5 \u2212 1 10 \u2212 3$ . | $ \Delta P L *$ . | 1 . |

Argon | 0.75–1.30 | GEMC | 65.428 | 26.382 | 5.950 | 0.373 | 0.447 | 15.437 |

CFC2F | 130.779 | 24.406 | 11.022 | 0.383 | 0.526 | 13.601 | ||

CFC1F | 95.114 | 7.227 | 2.585 | 0.267 | 0.381 | 10.330 | ||

Argon | 0.70–0.95 | GEMC | 196.474 | 0.426 | 1.261 | 0.031 | 0.530 | 10.077 |

CFC1F | 78.146 | 0.688 | 2.122 | 0.050 | 0.139 | 7.181 |

. | . | . | $ t Swap$ . | $ \Delta \rho V *$ . | . | $ \Delta P V *$ . | . | $ std ( n V )$ . |
---|---|---|---|---|---|---|---|---|

System . | $ T *$ . | Algorithm . | 1 . | $ \sigma 3 10 \u2212 3$ . | $ \Delta \rho L *$ . | $ \sigma 3 \u03f5 \u2212 1 10 \u2212 3$ . | $ \Delta P L *$ . | 1 . |

Argon | 0.75–1.30 | GEMC | 65.428 | 26.382 | 5.950 | 0.373 | 0.447 | 15.437 |

CFC2F | 130.779 | 24.406 | 11.022 | 0.383 | 0.526 | 13.601 | ||

CFC1F | 95.114 | 7.227 | 2.585 | 0.267 | 0.381 | 10.330 | ||

Argon | 0.70–0.95 | GEMC | 196.474 | 0.426 | 1.261 | 0.031 | 0.530 | 10.077 |

CFC1F | 78.146 | 0.688 | 2.122 | 0.050 | 0.139 | 7.181 |

. | . | $ t Swap$ . | $ \Delta \rho V *$ . | . | $ \Delta x V$ . | . | $ std ( n V )$ . |
---|---|---|---|---|---|---|---|

$ T *$ . | Algorithm . | 1 . | $ \sigma 3 10 \u2212 3$ . | $ \Delta \rho L *$ . | $ mol \u2009 mol \u2212 1 100$ . | $ \Delta x L$ . | 1 . |

0.75 | GEMC | 1,913.2 | 6.082 | 82.557 | 0.105 | 1.472 | 13.573 |

CFC1F | 50.2 | 6.911 | 52.737 | 0.234 | 2.167 | 18.403 | |

1.00 | GEMC | 43.3 | 4.274 | 101.496 | 0.938 | 1.854 | 18.073 |

CFC1F | 36.8 | 7.611 | 45.954 | 0.793 | 1.489 | 20.238 |

. | . | $ t Swap$ . | $ \Delta \rho V *$ . | . | $ \Delta x V$ . | . | $ std ( n V )$ . |
---|---|---|---|---|---|---|---|

$ T *$ . | Algorithm . | 1 . | $ \sigma 3 10 \u2212 3$ . | $ \Delta \rho L *$ . | $ mol \u2009 mol \u2212 1 100$ . | $ \Delta x L$ . | 1 . |

0.75 | GEMC | 1,913.2 | 6.082 | 82.557 | 0.105 | 1.472 | 13.573 |

CFC1F | 50.2 | 6.911 | 52.737 | 0.234 | 2.167 | 18.403 | |

1.00 | GEMC | 43.3 | 4.274 | 101.496 | 0.938 | 1.854 | 18.073 |

CFC1F | 36.8 | 7.611 | 45.954 | 0.793 | 1.489 | 20.238 |

Simulations with argon as a pure Lennard-Jones fluid were conducted for $ 0.75 < T * < 1.3$, specifically at 0.75, 1.00, 1.15, 1.25, and 1.30. The results of these simulations, which used initial conditions discussed in the literature, are summarized in Table II.^{10}

The coexisting densities and pressures are in good agreement with results from the Thol equation of state.^{37} The performance of the CFC1F algorithm is much better than the CFC2F algorithm, because the time needed for one successful swap is much less. Simulations for temperatures in the range $ 0.7 < T * < 0.95$^{24} were conducted for argon. For these lower temperatures, the time needed for one successful swap using the CFC1F algorithm is much less than the conventional GEMC algorithm.

### D. Asymmetric LJ mixture

The CFC1F algorithm was further investigated for three LJ mixtures at three pressures for each mixture. The mixture with the parameter ratios $ \sigma 22 / \sigma 11 = 0.5$ and $ \u03f5 22 / \u03f5 11 = 0.5$ is discussed here. The other two mixtures are discussed in Ref. 34. The results of the simulations are summarized in Table III and visualized in Fig. 8.

The conventional GEMC algorithm performed well at high temperatures ( $ T * = 1.00$) and needed a time equivalent to approximately 40 displacements for the asymmetric mixture to facilitate one successful swap. The performance of the algorithm clearly decreased at lower temperatures, and this effect is especially pronounced for asymmetric mixtures. The reason could be that the smaller component is two times smaller than the larger component. Thus, cavities suitable for an insertion of the larger component are easily filled with the smaller component, and the acceptance probability of the larger component decreases. Consequently, the computation time increased dramatically, and it took a time equivalent to approximately 2,000 attempted translations to obtain one successful swap.

In contrast to the conventional GEMC method, the CFC1F approach performed equally well for all investigated states. Although the GEMC algorithm performed better at high temperatures because of the higher computational effort of the CFC1F algorithm, it took a normalized time between 30 and 140 to facilitate one successful swap.

The t-PR equation of state^{38} predicted the molar fractions for those simulation points sufficiently well, and the initial guesses for densities yielded a reasonable computation time.

In brief, the GEMC and the CFC1F algorithms were in qualitative and quantitative agreement compared to the averaged deviation from the BWR equation of state. The deviation of the vapor density, $ \rho V *$, was on the order of $ 10 \u2212 3 \u2009 \sigma 3$, and the deviation of the liquid density, $ \rho L *$, was on the order of $ 10 \u2212 2 \u2009 \sigma 3$. The deviations of the molar fractions were of the order of 1%–5% for the simulated mixtures. The standard deviation of the number of molecules ranged between 10 and 30 for all algorithms and simulated states. This indicates that all the algorithms were similarly able to maintain stable box populations, which is required for efficient sampling.

## VII. DISCUSSION

We found that the CFC1F algorithm, which uses only one fractional, requires substantially less computational time, compared to CFC2F. The acceptance probability of insertions using the CFC1F algorithm ranged from 5% to 10%. However, more cycles are needed in some cases because CFC1F converges more slowly, compared to the traditional GEMC algorithm. A possible reason for the slower convergence could be that the average energy change due to a CFC move, compared to the GEMC swap move, is smaller. Therefore, uphill CFC moves are also more likely to be accepted, negatively affecting the speed at which the chemical potential equilibrates. In contrast, in the GEMC algorithm, the chemical potential and the energy contribution of a single molecule are directly related. The CFC1F algorithm involves less computational effort than the CFC2F algorithm, because the boxes are not coupled by the fractionals. Hence, CFC1F is a good choice for systems at low temperatures or asymmetric mixtures.

We conclude that it is better to use the conventional GEMC algorithm at moderate to high temperatures for pure LJ fluids. For low temperatures, dense phases, asymmetric mixtures, or complex molecules, the CFC1F algorithm is a good choice, because the speed of the expansion and contraction of the fractionals can be adjusted, ensuring reasonable acceptance probabilities.

## VIII. SUGGESTED PROBLEMS

Based on the results we have discussed and the simulations in Refs. 1, 10, 11, 12, 27, 38, and 40–42, we recommend the following simulation protocols.

The initial conditions should be chosen so that in the equilibrated system approximately 1/5 of the particles are in the vapor phase. Given that the molecules are moved between the boxes during the simulations, it is necessary to do more simulations if your initial guess is not accurate.^{40} For simulations of mixtures, the initial molar fractions should be equal for both boxes, $ x i , L = x i , V$.^{26} The values for the molar fractions should be chosen to be approximately in the center of the two-phase area in the vapor pressure (Pxy) diagram^{43} to assure enough particles in both phases.

The number of attempted displacements is initially set equal to the total number of molecules in the system, and one volume change is attempted per cycle. Although the number of swaps can be adjusted during the equilibration phase, recommended initial values are listed in Table IV.

$ T *$ | 0.75 | 1.00 | 1.15 | 1.25 | 1.30 |

$ n swaps / N tot$ | 1.50 | 0.20 | 0.10 | 0.04 | 0.04 |

$ T *$ | 0.75 | 1.00 | 1.15 | 1.25 | 1.30 |

$ n swaps / N tot$ | 1.50 | 0.20 | 0.10 | 0.04 | 0.04 |

Panagiotopoulos first proposed a sequence of displacements, volume changes, and swap moves.^{10} This order was later compared with a random choice of these trial moves.^{1} Because no significant differences in the simulation results were observed and the random arrangement of the trial moves ensures microscopic reversibility,^{11,12,44,45} random trial moves are recommended.

For the CFC2F algorithm, it is recommended that 50% of the attempted trial moves be displacements, 49% be *λ* moves, and 1% be volume changes. For the CFC1F algorithm, it is suggested that 49.5% of the trial moves be displacements, 49.5% be continuous fractional component trial moves, and 1% be volume changes.

We found that $ N tot / 4$ Widom insertions are sufficient because the energy contribution due to the insertion of an additional molecule is calculated each time a conventional GEMC swap is attempted. Because the alternative swap algorithms replace the GEMC swap, more Widom insertions are necessary if an alternative insertion algorithm is used. If CFC is used, fractionals are not considered in the energy calculation of Widom insertions. The range of the trial moves for the displacements, volume changes and the *λ* moves should be adjusted to yield an acceptance probability of $ \u2248 50 %$. For the displacements, $ 0.59 \u2009 \sigma $ is a reasonable initial value.

The initial maximum trial move for volume changes was set to 1% of the volume of the liquid box because the intermolecular forces in the liquid box are stronger than in the vapor box.

^{10}

^{27}In practice, we set the cut-off radius to half the box length, and it was dynamically updated if the box volume changed.

*Problem 1: Your first GEMC simulations*. We suggest starting with the GEMC algorithm, applied to argon as a pure Lennard-Jones fluid at $ T * = 1.0$. Simulations at $ T * = 1.0$ converge quickly, because swaps are sufficiently successful and $ T *$ is sufficiently far from the critical point. Choose $ \rho * = 0.06$, reduced liquid density $ \rho * = 0.60$, 150 molecules in the vapor box, and 350 molecules in the liquid box.^{10} Use 2,000 MC cycles to melt the initial face centered cubic grid, 2,000 MC cycles to equilibrate the system, and 8,000 MC cycles to obtain a sufficient number of states for the ensemble averages. Check the plausibility of your results by analyzing a Frenkel plot, as illustrated in Fig. 9.

As shown in Fig. 9, the successful simulation quickly converges to its equilibrium configuration and samples the phase space efficiently. The equilibrium configuration is sampled much more frequently (see the probability indicator in the bar on the right-hand side of the plot for each phase). In contrast, the unsuccessful simulation resulted in a state where most molecules populate only one of the boxes.

After this initial check, assess the time series of the number of molecules per box, the box volumes, the box pressures, and the chemical potential. If the running average of any of these quantities is not well defined, the duration of the simulation needs to be extended. You should find that the densities are $ \rho V * = 0.03$ and $ \rho L * = 0.70$ and that the phase pressure is $ P * \u2248 0.02$. The initial conditions of the unsuccessful simulation led to a configuration where one of the boxes became almost empty. One way to find the tipping point between successful and unsuccessful initial conditions is to keep the initial densities constant and vary the number of molecules per box. For example, we can place 480 molecules in the vapor box and only 20 molecules in the liquid box. Another option is to keep the number of molecules per box constant and vary the initial molar densities. Use one of these suggestions or your own choices and conduct short simulations (2,000 warm-up cycles, 2,000 equilibration cycles, 500 production cycles) to find that tipping point and evaluate the recommendations in Sec. VIII. If you use the initial conditions suggested in Ref. 10 for $ T * = 1.0$, you will find that the system starts to become unstable when less than 60 molecules are placed in the vapor box.

*Problem 2: Using the continuous fractional component algorithm*. Below a certain temperature you will find that most of the time is consumed by unsuccessful attempts to swap molecules between the two phases, and the continuous fractional component algorithm becomes useful. What is the value of $ T *$ at which it makes sense to use the continuous fractional component algorithm? Compare and analyze the course of simulations conducted with the GEMC and the CFC1F algorithm. What are the most significant differences for the results you obtain? Example results are displayed in Fig. 10, which shows the results of the conventional GEMC algorithm, the CFC1F algorithm using the same overlap as the GEMC algorithm, and the CFC1F algorithm with an overlap of 1 Å. As expected, the acceptance probability decreases with decreasing temperature. This trend is much more pronounced for the GEMC algorithm than for the CFC1F algorithm. Apart from that, the acceptance probability is significantly higher when the overlap is set to 1 Å. If the overlap is chosen to be too large, moves of a fractional with a low coupling value of *λ* would be rejected by the overlap check, while a Metropolis check would have accepted them.

Does it make sense to further decrease the overlap radius for the continuous fractional component algorithm or might there be a limiting value below which it does not make sense to further decrease it? Which other criterion could be used to determine an overlap radius besides the interaction energy?

*Problem 3: Simulating mixtures*. To simulate mixtures, the initial conditions need to be chosen with even more care. In addition to the total system volume, the phase composition also needs to be chosen so that the mixing point is within the two-phase area; otherwise, the lever rule^{43} will leave a phase with little or no molecules. Use an appropriate equation of state (e.g., t-PR^{38}) to estimate the initial molar fractions. Choose the CFC1F algorithm if the acceptance probability is very low, for example, when size-asymmetric mixtures or low temperatures are simulated. The acceptance probability of the larger component might be much different than that of the smaller component. For example, the larger molecules might create cavities which can accommodate insertions of the smaller molecules. A larger molecule is also more likely to overlap with any other molecule in the system. Consequently, the chemical potential might not be properly equilibrated because a sufficient number of accepted swaps is required. These results are seen in Fig. 11, where the acceptance probability of the larger component and the smaller component of each binary mixture is compared for the GEMC and CFC1F algorithm at $ T * = 1.0$.

For the CFC1F algorithm, an acceptance probability for a swap move on the order of 5% to 10% was observed. This higher acceptance probability is especially pronounced for the most asymmetric mixture due to the fact that the maximum change of the coupling factor, *λ*, which controls the growth and shrinking of the fractional, is adjusted separately for each component, whereby a reasonable acceptance probability for both components, even at significant size differences, is obtained. When deciding whether it is feasible to use the CFC1F algorithm, the increased computational time needs to be compensated. Example results for the normalized time can be found in Table III. You will find that it is reasonable to use the CFC1F algorithm at low temperatures such as $ T * = 0.75$, while the GEMC algorithm performs better at moderate to high temperatures such as $ T * = 1.0$. For mixtures, it is reasonable to use the CFC1F algorithm for low temperatures $ T * < 0.9$, or for asymmetric mixtures, independent of the simulated temperature.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

## REFERENCES

*Computer Simulation of Liquids*

*CRC Handbook of Chemistry and Physics*

*Applications of Molecular Simulation in the Oil and Gas Industry: Monte Carlo Methods*

*Atkins' Physical Chemistry*

*Understanding Molecular Simulation*