Molecular exchange Monte Carlo: A generalized method for identity exchanges in grand canonical Monte Carlo simulations

In Monte Carlo simulations in the grand canonical ensemble (GCMC), the chemical potential, volume, and temperature are fixed (μVT = constant). Sampling of phase space is achieved through a variety of trial moves, such as displacement and molecule insertion and deletion. For complex molecular typologies, additional trial moves, such as rigid body rotation and configurational-bias regrowth,^1,2 may be included to improve the sampling of conformational degrees of freedom. During the course of the simulation, the conjugate variables N (number of molecules) and E (potential energy) fluctuate. Because GCMC allows for the simulation of an open system, it has been used extensively to study the adsorption of gases in porous materials.^3–6 When combined with histogram-reweighting methods,^7,8 GCMC simulations provide precise predictions of vapor-liquid equilibria for pure fluids and mixtures^9,10 and have been used to determine critical micelle concentrations for model surfactants.¹¹ 

Perhaps the greatest challenge with GCMC simulations, however, is achieving a sufficient number of accepted molecular insertion/deletion moves to ensure adequate sampling of phase space. Therefore, significant effort has been expended to develop algorithms that improve the acceptance rate for molecule insertions and deletions. Biasing methods, such as rotational, energy, and cavity-bias, were used to improve the efficiency of simulations for the adsorption of benzene and p-xylene in silicalite.¹² The introduction of configurational-bias Monte Carlo (CBMC) enabled the successful simulation of chain molecule adsorption in zeolites,¹³ which was followed by the coupled-decoupled¹⁴ and reservoir methods,^15,16 which extended the complexity of systems that could be simulated to include molecules with branch points and rings.

These aforementioned biasing methods have greatly extended the complexity of systems that may be simulated with GCMC simulations; however, at high densities and low temperatures, the acceptance rate for molecule transfers is still unacceptably low due to the difficulty in finding a favorably sized cavity to insert a molecule. For example, in simulations of branched alkanes, acceptance rates for molecule transfers at 0.7T_c were approximately 0.3%.⁹ Others have sought to address these issues through the use of cavity-bias,^17–19 to identify favorable locations to attempt molecule insertions, or continuous fractional component Monte Carlo (CFCMC),^20,21 and expanded ensembles,^22,23 where molecules are gradually inserted while the system is allowed to relax locally to minimize steric and energetic penalties due to molecule insertion.

For mixtures, a straight forward approach is to introduce a trial move where the identity of one molecule is changed to that of another.²⁴ The benefit of such a move is that steric overlaps are reduced significantly, leading to enhanced acceptance for the particle exchange. The identity exchange move has been used in many simulations of single particles in various ensembles, such as semi-grand,^25,26 Gibbs,^24,27,28 and grand canonical.^29–31 The methodology has been extended to allow for the exchange of multiple solvent molecules with a polymer chain composed of solvent monomers without changing the coordinates of either a polymer or solvent.³¹ For the simulation of mixtures of colloids and solvent, it is necessary to swap a large colloid particle for multiple smaller solvent particles. By swapping multiple solvent particles, it is possible to create large enough voids such that a reasonable acceptance rate may be obtained for the insertion of colloid particles.^29,30 For the exchange of a large particle with multiple small ones, Vink et al. used simple random insertions to determine the coordinates for the solvent particles. When inserting a large number of solvent particles, the potential for overlap increases, reducing the efficiency of the method. To address this issue, Kindt introduced the idea of “solvent repacking” for a two-dimensional hard-disk and size asymmetric three-dimensional Lennard-Jones systems, where configurational-bias was used to determine the positions of solvent particles in the large-small particle identity exchange.^29,32 The identity exchange move has also been used in Gibbs ensemble Monte Carlo simulations for molecular systems.^33–36 

The previously described methods for identity exchange were generally applicable to only the special cases for which they were developed, e.g., single particle exchanges,²⁶ a polymer composed of solvent monomers,³¹ or large hard particles or disks in a solvent of smaller hard particles.^29,30 These methods are difficult to generalize to molecular systems of arbitrary molecular topology, and their computational performance is expected to be highly correlated with the type of system for which the move was originally developed. To address these issues, a generalized identity exchange move for simulations in the grand canonical ensemble, referred to as Molecular Exchange Monte Carlo (MEMC), is presented which works for systems of any molecular topology. Three different approaches for the insertion of the large molecule are presented. A derivation of acceptance criteria and the algorithms for performing the MEMC move are provided in Sec. II for each of the three approaches, while the detailed computational procedure and mathematical calculations are included in the supplementary material. The utility of the three methods and their computational efficiency is illustrated for selected binary mixtures in Sec. IV. The key findings of the work are summarized in the Conclusions.

To describe the MEMC move in the grand canonical ensemble, it is helpful to consider the case of a large molecule that is exchanged with multiple smaller molecules. However, the methods may be applied without modification to the exchange of molecules of similar size. The original state is called the old state, while the state created by the attempted exchange move is called the new state. For a given configuration, with N_L large and N_S small molecules, an “insertion move” is an attempt to exchange one large molecule with N_EX small molecules inside a predefined exchange sub-volume V_EX, and a “deletion move” is an attempt to exchange N_EX small molecules for a large one. The exchange sub-volume is defined as an orthogonal box, where the length of the box in the x and y dimensions is set to the same values for simplicity and the z dimension is set independently. If desired, all three sub-volume box dimensions could be set independently. An orthogonal sub-volume is used instead of a cube or sphere to accommodate large molecules with different aspect ratios. Depending on the method used, the orientation of the exchange sub-volume z-axis may also be varied. Although not used in this work, it is also possible to optimize N_EX and V_EX “on the fly” during a simulation to maximize the acceptance rate.

The acceptance criterion for a molecular exchange move that satisfies the detailed balance equation is written as

where K(i → j) is the flux of probability from state i to state j. The probability flux is equal to the product of the probability of finding the system in state i, the probability of generating a move that takes state i to state j, and the probability of accepting the move

Based on the detailed balance Eq. (1), the ratio of the probability of accepting the move from old → new to that of its reverse move new → old is

In the deletion move, where one large molecule is exchanged for N_EX small molecules, the ratio of the probability of being in the configuration new to the probability of being in the configuration old is

where $β=1/kBT$ and $μL$ and $μS$ are the imposed chemical potentials of large and small molecules, respectively. U(old) and U(new) are the potential energies of the system in configuration old and configuration new, respectively.

For the insertion move, where N_EX small molecules are exchanged for one large molecule, the ratio of the probability of being in the configuration new to the probability of being in the configuration old is

The probability of generating the new state, for both insertion and deletion of the large molecule, is given by the product of the probability of locating the center of the exchange sub-volume at a particular point within the simulation box, the probability of choosing N_EX particular small molecules, the probability of choosing a particular large molecule, the probability of generating trial configurations for N_EX small molecules, and the probability of generating trial configurations for the large molecule,

Depending on how the center of the exchange sub-volume is located, the molecules to be exchanged are chosen, and depending on how trial positions are generated, different algorithms to perform the MEMC move may be devised.

For the large molecule insertion move, the exchange sub-volume V_EX with a random geometric center and a random orientation is defined within the simulation box. For a large molecule deletion move, the geometric center of V_EX is located at the center of mass (COM) of the selected large molecule and its z-axis is aligned with the backbone of the large molecule. See Fig. 1 for more details.

The algorithm for the insertion of a large molecule after the deletion of small molecule(s) is as follows:

1. Define an orthogonal exchange sub-volume V_EX, with its geometric center located randomly within the simulation box of volume V (with the probability proportional to V⁻¹) and a random orientation. Determine the total number of small molecules within the exchange sub-volume (N_S,VEX) based on their center of mass.

2. Reject the move if $NS,VEX<NEX$⁠, otherwise continue.

3. Select N_EX small molecules out of N_S,VEX found in the exchange sub-volume with the probability of $NEX!(NS,VEX−NEX)!/NS,VEX!$⁠.

4. Repeat the steps a and b for N_EX cycles (i = 1, 2,…, N_EX) to delete the selected small molecules.

   • Generate j − 1 random trial positions for the COM of the ith small molecule within the exchange sub-volume V_EX. The original position of the COM of the ith small molecule will be included as the jth term.

   • For each trial COM position p, generate k random trial orientations around the molecule’s COM (except the jth COM, where k − 1 random trial orientations are generated and the original orientation of the molecule will be included as the kth term) and calculate the Rosenbluth weight $Wi,old=∑p=1j∑r=1k⁡exp−βUi,p,r$⁠, where U_i,p,r is the interaction energy of the ith molecule to be removed in position p and orientation r with all other molecules, excluding those removed in the earlier cycles of the move. Finally, remove the molecule from the simulation box. Calculate $Pi,old=exp−βUi,j,kWi,old,$ where U_i,j,k is the interaction energy of the ith small molecule at its original COM position and orientation with all other molecules remaining in the simulation box. P_i,old is the probability of inserting the ith small molecule back in its original configuration in the reverse move (new → old).

5. Insert the COM of the large molecule at the geometric center of the exchange sub-volume V_EX and align the backbone of the large molecule with the z-axis of the exchange sub-volume. Generate k random trail orientations for the large molecule around the z-axis (two-dimensional rotation). Calculate the Rosenbluth weight $Wnew=∑r=1k⁡exp−βUr$⁠, where U_r is the interaction energy of the inserted large molecule at orientation r with all other molecules in the simulation box.

6. Select one of the generated trial configurations with the probability $Pnew=exp−βUrWnew$ and insert the large molecule.

The algorithm for the deletion of a large molecule and the subsequent insertion of small molecule(s) is as follows:

1. Select a large molecule out of N_L large molecules within the simulation box with the probability of 1/N_L.

2. Define an orthogonal exchange sub-volume with its geometric center placed at the COM of the selected large molecule and its z-axis aligned with the backbone of the large molecule. Determine the number of small molecules N_S,VEX within the exchange sub-volume.

3. Generate k − 1 random trial orientations for the large molecule around the z-axis of the sub-volume. The original orientation will be included as the kth term in the Rosenbluth weight. The Rosenbluth weight is calculated as $Wold=∑r=1k⁡exp−βUr$⁠, where U_r is the interaction energy of the large molecule in orientation r with all other molecules in the simulation box. Calculate the probability $Pold=exp−βUkWold$⁠, where U_k is the interaction energy of the large molecule at the original orientation with all other molecules in the simulation box. P_old is the probability of inserting the large molecule at its original configuration in the reverse move (new → old). Then remove the large molecule from the simulation box.

4. Repeat the steps a → c for N_EX cycles (i = 1, 2,…, N_EX) to insert the small molecules with the probability of $NEX!/VEXNEX$⁠.

   • Generate j random trial positions for the COM of the ith small molecule within V_EX.

   • For each trial position p, generate k random trial orientations around the molecule’s COM (three-dimensional rotation) and calculate the Rosenbluth weight $Wi,new=∑p=1j∑r=1k⁡exp−βUi,p,r$⁠, where U_i,p,r is the interaction energy of the ith inserted small molecule at position p and orientation r with all the other molecules, including those added in the earlier cycles of the move.

   • Pick one of the generated trial configurations with the probability $Pi,new=exp−βUi,p,rWi,new$ and insert the small molecule.

Based on the two algorithms described above, for the large molecule insertion, the ratio of the probability of generating the move new (N_L + 1, N_S − N_EX) → old (N_L, N_S) to that of the reverse move is

Simplifying Eq. (7) and substituting into Eq. (3) produces the acceptance criteria for the large molecule insertion,

For the large molecule deletion move, the ratio of the probability of generating the move new (N_L − 1, N_S + N_EX) → old (N_L, N_S) to that of the reverse move is

Simplifying Eq. (9) and substituting into Eq. (3) produces the acceptance criteria for the large molecule deletion move

The energy difference between configurations new and old, U(new) − U(old), does not appear directly in the acceptance criteria because their Boltzmann weight is already included in the probabilities used for selecting the position of the molecules.

The acceptance criterion derived for ME-1 is identical to the one introduced by Vink and Horbach.³⁰ This move performs well for binary mixtures with low concentrations of large molecules. However, the acceptance rate of the move decreases significantly as the concentration of the large molecules increases, and the chance of finding N_EX small molecules in the exchange sub-volume becomes very low. To address this limitation, ME-2 was developed.

In ME-1, for the insertion of a large molecule, the exchange sub-volume V_EX is defined with a random orientation and position. However, as the mole fraction of small molecules decreases, the required number of small molecules is frequently not available within the exchange sub-volume. Therefore, a large fraction of the attempted ME-1 moves will be rejected. In the ME-2 approach, the geometric center of V_EX is placed on the COM of a randomly selected small molecule. If the small molecule is monoatomic, the orientation of V_EX is assigned randomly, otherwise its z-axis is aligned with the backbone of the small molecule. The large molecule deletion is identical to ME-1. An illustration of the ME-2 algorithm is provided in Fig. 2.

The algorithm for the insertion of a large molecule after the deletion of small molecule(s) is as follows:

1. Select one molecule out of N_S small molecules in the simulation box with the probability of 1/N_S. This molecule will be the last molecule to be removed from the system.

2. Define V_EX with its geometric center placed at the COM of the small molecule selected in step 1. The z-axis of the exchange sub-volume is aligned with the backbone of the small molecule. If the small molecule is monoatomic, the orientation of V_EX is assigned randomly. Determine the number of small molecules N_S,VEX within V_EX (N_S,VEX includes the molecules selected in step 1).

3. Reject the move if $NS,VEX<NEX$⁠, otherwise continue.

4. Select N_EX − 1 small molecules out of N_S,VEX − 1, with the probability $(NEX−1)!(NS,VEX−NEX)!/(NS,VEX−1)!$⁠.

5. Repeat steps a and b of the large molecule insertion move of ME-1 for N_EX − 1 cycles (i = 1, 2,…, N_EX − 1) to delete the selected small molecules.

6. For the last small molecule to be deleted, generate k − 1 random trial orientations around the z-axis of the sub-volume. The original orientation will be included as the kth term in the Rosenbluth weight. The Rosenbluth weight is calculated from $WNEX,old=∑r=1k⁡exp−βUNEX,r$⁠, where U_NEX,r is the interaction energy of the last small molecule in orientation r with all other molecules in the simulation box. Finally, remove the last small molecule from the simulation box and calculate $PNEX,old=exp−βUNEX,kWNEX,old,$ where U_NEX,k is the interaction energy of the last small molecule at its original configuration with all other molecules remaining in the simulation box. P_i,old is the probability of inserting the ith small molecule back in its original configuration in the reverse move (new → old).

7. Insert the large molecule according to steps 5 and 6 of ME-1.

The algorithm for the deletion of a large molecule and subsequent insertion of small molecule(s) is as follows:

1. Follow steps 1-4 of the ME-1 large molecule deletion move.

2. Insert the COM of the first small molecule at the geometric center of V_EX and align its backbone with the z-axis of the exchange sub-volume. Generate k random trial orientations around the z-axis of the sub-volume. Calculate the Rosenbluth weight $W1,new=∑r=1k⁡exp−βU1,r$⁠, where U_1,r is the interaction energy of the first small molecule inserted at orientation r with all other molecules in the simulation box.

3. Select one of the trial orientations with the probability $P1,new=exp−βU1,rW1,new$⁠.

4. Repeat steps a → c of the large molecule deletion move of ME-1 for N_EX − 1 cycles (i = 2,…, N_EX) to insert the small molecules with the probability $NEX−1!/VEXNEX−1$⁠.

Based on the two algorithms described above, for the large molecule insertion move, the ratio of the probability of generating the move new (N_L+1, N_S − N_EX) → old (N_L, N_S) to that of the reverse move is

Simplifying Eq. (11) and substituting into Eq. (3) results in the acceptance criterion for the large molecule insertion move

For the large molecule deletion move, the ratio of the probability of generating configuration $newNL−1,NS+NEX→oldNL,NS$ to that of the reverse move is

Simplifying Eq. (13) and substituting into Eq. (3) results in the acceptance criterion for the large molecule deletion move,

If N_EX = 1, the acceptance criteria given in Eqs. (13) and (14) simplify to that of the standard identity-exchange acceptance move,²⁶ 

For the large molecule insertion move in ME-2, the large molecule is inserted as a rigid body and its backbone is aligned with the z-axis of the V_EX. This move performs well for large molecules with a straight backbone. However, the acceptance rate decreases for a large molecule with nonlinear geometry as it becomes significantly more difficult to fit a complex rigid body into the void space created after deleting the small molecule(s). Therefore, a modification to ME-2 was developed to address this limitation.

In the ME-3 algorithm, a predefined atom of the large molecule is first placed at the geometric center of V_EX and the molecule is built segment by segment using the coupled-decoupled configurational-bias Monte Carlo (CBMC) algorithm.¹⁴ For the large molecule deletion move, the exchange sub-volume is defined with a random orientation, with its geometric center placed at the same predefined atom of the large molecule to be deleted. Next, the Rosenbluth weight W_old of the large molecule is calculated. Insertion and deletion of N_EX small molecules is identical to the ME-2 method. Figure 3 illustrates the ME-3 algorithm.The algorithm for the insertion of a large molecule after deletion of small molecule(s) is as follows:

1. Select one molecule out of N_S small molecules in the simulation box with the probability 1/N_S. This molecule will be the last molecule to be removed from the system.

2. Define an orthogonal exchange sub-volume V_EX with a random orientation and its geometric center placed at the COM of the small molecule selected above. Then determine the number of small molecules N_S,VEX within V_EX (N_S,VEX includes the molecule selected in step 1).

3. Repeat steps 3-6 of the ME-2 method to delete N_EX small molecules from the simulation box.

4. Insert the predefined atom of the large molecule at the center of V_EX and perform coupled-decoupled configurational-bias Monte Carlo to grow the large molecule segment by segment. Calculate the Rosenbluth weight W_new.

5. Insert the large molecule by selecting one of the generated trial configurations with the probability P_new.

The algorithm for the deletion of a large molecule and subsequent insertion of small molecule(s) is as follows:

1. Within the simulation box of volume V, pick one large molecule out of N_L with the probability of 1/N_L.

2. Define an orthogonal exchange sub-volume V_EX with a random orientation and place its geometric center at the predefined atom of the selected large molecule. Determine the number of small molecules N_S,VEX within the exchange sub-volume.

3. Perform coupled-decoupled CBMC for the large molecule and calculate the Rosenbluth weight W_old and P_old.

4. Repeat steps 2-4 of ME-2 to insert N_EX small molecules within V_EX.

The forward to reverse probability ratios for generating the large molecule insertion and the large molecule deletion moves are identical to those given in Eqs. (11) and (13), respectively. The acceptance criteria for the ME-3 algorithm are identical to those of ME-2 and are given by Eqs. (12) and (14).

The three molecular exchange algorithms described in Sec. II were implemented in the development version of GPU Optimized Monte Carlo (GOMC), which is available to the public on GitHub.³⁷ GOMC is an object-oriented Monte Carlo simulation engine, capable of performing simulations in canonical, isobaric-isothermal, and grand canonical ensembles, as well as Gibbs ensemble Monte Carlo. GOMC is designed for the simulation of complex molecular topologies and supports a variety of potential functions, such as Lennard-Jones and Mie potentials. Coulomb interactions are also supported via the Ewald summation method.³⁸ GOMC is capable of parallel computation, on either multicore CPUs or GPUs.

Phase diagrams were determined from histogram-reweighting Monte Carlo simulations in the grand-canonical ensemble.³⁷ A cubic box size of 25 Å × 25 Å × 25 Å was used for methane+ethane, methane+propane, methane+n-butane, and water+impurity. For perfluorobutane+n-butane and methane+n-pentane, a box size of 30 Å × 30 Å × 30 Å was used, while for 2,2,4-trimethylpentane+neopentane a box size of 40 Å × 40 Å × 40 Å was used. Initial configurations were generated with Packmol.³⁹ Psfgen was used to generate coordinate (^*.pdb) and connectivity (^*.psf) files.⁴⁰ Potentials were truncated at 10 Å, and analytical tail corrections were applied.⁴¹ To enhance the acceptance rate for molecule insertions, the coupled-decoupled configurational-bias Monte Carlo (CBMC) algorithm was used.¹⁴ For all liquid phase simulations, unless otherwise noted in Sec. IV, configurational-bias parameters were 100 angle trials, 100 dihedral trials, 10 trial locations for the first site, and 8 trial locations for secondary sites. For standard GCMC simulations, a move ratio of 20% displacements, 10% rotations, 10% regrowth, and 60% molecule transfers was used. For simulations that included the molecular exchange move, 30% molecular exchanges were performed with a corresponding reduction in the percentage of attempted molecule transfers.

Uncertainties used in the calculation of the statistical efficiency of the methods were calculated as the standard deviation determined from five unique simulation trajectories, each started from a unique initial configuration and random number seed. All simulations, except those used to generate phase diagrams, were run for 2 × 10⁷ Monte Carlo steps (MCS), without equilibration period. Simulations used to generate phase diagrams were run for 5 × 10⁷ MCS with a 5 × 10⁶ MCS equilibration period. Every 200-500 MCS, the instantaneous state of the system (N₁, N₂, E) was saved as a histogram. Every one million MCS, the natural log of distribution of large particles ln(P_N) for each simulation was determined, and the standard deviation and efficiency were calculated for each binary system for a variety of compositions along the vapor-liquid coexistence curve. Calculations were performed on one core of an Intel Xeon E5-4627v4 2.6 GhZ CPU.

The efficiency was calculated using the calculated standard deviation and the CPU time,

where $σ$ is average uncertainty in the natural log of large particle distribution and s is the CPU time in seconds.

In this section, a number of examples are provided to illustrate the effect of molecular exchange moves on the statistical sampling in grand canonical histogram reweighting Monte Carlo simulations. Mixtures simulated include perfluorobutane+n-butane and methane+ethane, +propane, +n-butane, and +n-pentane. Additional calculations were performed to generate pure fluid phase diagrams for water and 2,2,4-trimethylpentane to demonstrate the utility of the method and to provide comparisons to prior work.^36,42,43 For binary mixture phase diagrams, all calculations were performed at temperatures below 0.7T_c. For pure fluid phase diagrams, calculations were performed from the critical temperature to 0.44T_c – 0.51T_c. Performing grand canonical Monte Carlo simulations, using standard configurational-bias methods,¹⁴ below 0.7T_c is a challenging task, and therefore a good test to evaluate the improvement in sampling of phase space provided by the proposed algorithms.

Methane+n-alkane systems are well studied and extensive experimental data may be found in the literature.^44–51 In general, the determination of vapor-liquid coexistence for these systems at temperatures above 0.7T_c can be done using standard configurational-bias methods in grand canonical or Gibbs ensemble Monte Carlo simulations.^9,52–54 However, below 0.7T_c, acceptance rates for the insertion of n-alkanes into a liquid phase drop to approximately 0.1%, which necessitates long simulations to obtain convergence of the simulations. In this section, the effect of the three ME algorithms on the convergence of grand canonical Monte Carlo simulations is assessed for mixtures of methane+ethane, +propane, +n-butane, and +n-pentane, and the effectiveness of performing a two for one exchange is evaluated.

The methane+n-butane mixture is presented first as an example of the validation process used in the development of the molecular exchange methods. Grand canonical Monte Carlo (GCMC) simulations were performed for a variety of temperatures, chemical potentials, and move ratios using both standard configurational-bias insertions/deletions and the ME-1, ME-2, and ME-3 methods. Probability distributions of states sampled during the simulation were collected and compared to reference distributions determined using standard configuration-bias insertions. An example of this is shown in Fig. 4 for gas (μ_butane = −2960, μ_methane = −2000) and liquid (μ_butane = −2840, μ_methane = −2000) phase simulations at 277 K. As expected, the probability distributions produced by the ME-3 algorithm are an exact match to the reference distributions. Additional data for the ME-1 and ME-2 algorithms are presented in Figs. S1 and S2 of the supplementary material.

In Fig. 5, the pressure vs. composition diagram for methane+n-butane at 277 K, predicted using both the coupled-decoupled configurational-bias method¹⁴ and the ME-3 algorithm, is shown. Interactions between molecules were described with optimized Mie potentials for phase equilibria.⁵² In addition to showing excellent agreement with experimental data,⁵¹ the ME-3 algorithm produced results that were statistically indistinguishable from standard configurational-bias insertions, providing further validation of the method.

In Table I, the acceptance rate for molecule transfers as a function of composition is presented for each methane+n-alkane binary mixture. Calculations were performed for liquid phase simulations along the coexistence curve at 186 K (methane+ethane), 213 K (methane+propane), 225 K (methane+n-butane), and 273 K (methane+n-pentane). The systems exhibit similar general trends, with acceptance rates climbing as the critical point of the mixture is reached. For $xmethane<0.5$⁠, acceptances rates for the insertion of the larger n-alkane using configurational-bias were less than 1%. When performing a one-to-one exchange, ME-3 was found to produce the largest improvement in acceptance rates for the large molecule, producing improvements of 2× for methane+n-pentane at x_methane = 0.7 to 70× for methane+ethane at x_methane = 0.1. The ME-2 algorithm also produced a significant enhancement in the acceptance rates for the insertion of the longer n-alkane, while the ME-1 algorithm was found to yield significantly lower acceptance rates than traditional configurational-bias insertions. Because the ME-2 algorithm uses a rigid swap and the COM of the large molecule is placed at the geometric center of the exchange sub-volume, only a fraction of the sub-volume is guaranteed to be empty. In most of the ME-2 exchanges, it is likely that some atoms from the large molecule will have overlaps with existing molecules, lowering acceptance rates compared to ME-3. The ME-3 algorithm uses the same initial placement for the central atom as ME-2, but grows the rest of the large molecule, allowing it to find more energetically favorable configurations than are possible through a rigid molecule insertion, leading to greater acceptance rates for the exchange move. As expected, the more similar the large and small molecules are in terms of excluded volume, the greater the success of the molecular exchange. It is also interesting to note that even for the highly asymmetric mixture of methane+n-pentane, acceptance rates for molecule transfers were improved substantially through the inclusion of the molecular exchange move.

The molecular exchange algorithm allows for trial moves where any number of small molecules may be exchanged for one large molecule. An example of this is shown in Table I, where acceptance rates are presented for the exchange of two methanes with one n-butane or n-pentane (N_EX = 2). For the ME-3 algorithm, acceptance rates are always lower than the one for one exchange, although this difference decreases as the chain length of the large molecule increases. Part of the decrease in the acceptance rate stems from the reduced probability of finding two methane molecules in the sub-volume to exchange at low methane concentrations. For ME-2, acceptance rates are slightly lower for the exchange of two methanes with one n-butane, compared to the one-to-one exchange. However, for the exchange of two methanes with one n-pentane, slight improvements in the acceptance rates were observed. The ME-1 algorithm shows a slight improvement in acceptance rates for the exchange of two methanes with one n-butane or n-pentane, although in all cases, acceptance rates for the ME-1 algorithm are lower than for configurational-bias insertions.

While the size of the sub-volume does not have an effect on the acceptance rates for the ME-2 and ME-3 algorithms for a one-to-one exchange, it was found to have an effect on the acceptance rates for the two-to-one exchange, as shown in Table I. Increasing the size of the sub-volume increases the probability that a second small molecule will be found within the sub-volume, leading to an increased overall acceptance rate for the MEMC move. Therefore, it is possible to optimize acceptance rates for the two-to-one exchange ratio by performing a series of short simulations for a range of sub-volume box lengths and by using a heuristic that the sub-volume should be large enough to contain the entire large molecule. For methane+n-butane, the optimum exchange sub-volume size for a two-to-one exchange was found to be 8.8 Å × 8.8 Å × 11.8 Å for ME-3 and ME-2, while for ME-1 it was 5 Å × 5 Å × 8 Å.

A more detailed analysis of the statistical uncertainty and efficiency for an exchange ratio of one-to-one is provided in Fig. 6 for the methane+n-butane mixture. A direct comparison between the efficiencies obtained for the one-to-one and one-to-two exchange ratios is presented in Fig. S3 of the supplementary material. Uncertainties were determined from probability distributions collected from liquid phase grand canonical Monte Carlo simulations performed along the vapor-liquid coexistence curve. For all mole fractions investigated, the ME-3 algorithm shows the fastest convergence of the n-particle probability distribution, converging in approximately half the number of Monte Carlo steps of ME-2. Both the ME-3 and ME-2 produce similarly converged probability distributions after 2 × 10⁷ MCS, with average uncertainties of approximately 0.05. The ME-1 algorithm and configurational-bias insertions show similar convergence properties. However, with 2 × 10⁷ MCS, each produced uncertainties that were approximately double than those of the ME-3 and ME-2 methods.

In Fig. 7, the probability distributions resulting from GCMC simulations with the various ME methods using an exchange ratio of one to one are presented for x_methane = 0.3, while data for other mole fractions are given in Figs. S4 and S5 of the supplementary material. The probability distributions resulting from GCMC simulations with the various ME methods using an exchange ratio of one to two are presented in Figs. S6–S8 of the supplementary material for a range of mole fractions. All MEMC methods converge to the same distribution. ME-3 shows rapid convergence, and within only 5 × 10⁶ MCS, the correct distribution is obtained. The ME-2 algorithm shows slightly slower convergence compared to ME-3 but is still more efficient than ME-1 or configurational-bias trial insertions.

The perfluorobutane+n-butane system is an interesting case study because of its large deviations from Raoult’s law, despite the fact that perfluorobutane and n-butane have very similar normal boiling points (270.96 K for C₄F₁₀ and 272.61 K for C₄H₁₀) and both are non-polar with similar molecular geometries. This system has been modeled in the past with Statistical Associating Fluid Theory Variable Range (SAFT-VR),⁵⁵ Perturbed Chain Statistical Associating Fluid Theory (PC-SAFT),⁵⁶ and Group Contribution Statistical Associating Fluid Theory Variable Range (GC-SAFT-VR),⁵⁷ which showed close agreement with experimental data.⁵⁸ Gibbs ensemble Monte Carlo simulations using an identity exchange move have been used to study liquid-liquid equilibria for n-heptane+perfluoheptane,⁵⁹ otherwise, grand canonical and Gibbs ensemble methods have rarely been applied to these kinds of mixtures. This is due, in part, to the difficulty in achieving an adequate number of accepted molecule transfers. For example, at 260 K, acceptance rates for the insertion of perfluorobutane in the neat liquid phase were approximately 0.075%.

In Fig. 8, the pressure vs. composition diagram for perfluorobutane+n-butane at 260 K, predicted using the ME-3 algorithm and the Mie potentials developed by our group,⁵² is shown. The force field for perfluorobutane was modified slightly from the original work to use a more accurate seven term cosine series, which is described in detail in the supplementary material. Using standard Lorentz-Berthelot combining rules^60,61 and no adjustable parameters for the cross interaction, very good agreement was achieved with experiment. The largest deviation results from the limitation in the united-atom force field for perfluorobutane, which over-predicts the vapor pressure at 260 K by approximately 0.1 bar.

To evaluate the effectiveness of the molecular exchange move with a one-to-one exchange ratio and an exchange sub-volume of 6 Å × 6 Å × 9 Å, acceptance rates, uncertainties in the probability distributions, and efficiencies produced from the grand canonical Monte Carlo simulations were determined for liquid phase simulations at selected points along the coexistence curve. The effect of various simulation parameters on the performance of the CBMC and MEMC acceptance rates and efficiencies was also evaluated for liquid phase simulations containing 50 mol. % n-butane and is shown in Fig. S9 of the supplementary material. Using the coupled-decoupled configurational-bias method,¹⁴ the probability of successfully inserting one perfluorobutane into a simulation box containing 10 mol. %, 50 mol. %, and 90 mol. % of n-butane was 0.073%, 0.026%, and 0.011%, respectively. The ME-1 algorithm increased acceptance rates approximately 4 times that of standard trial insertions for $xbutane>0.50$⁠; however, for lower concentrations of n-butane, no improvement was observed. For the ME-2 algorithm, acceptance rates of 4.92%, 4.17%, and 3.15% were obtained, while for ME-3, acceptance rates were 3.52%, 2.73%, and 1.69%, respectively. For this system, the ME-2 algorithm produces the best acceptance rates because it works by aligning the backbone of perfluorobutane with the cavity left by the leaving n-butane. Acceptance rates were slightly lower for ME-3 since it grows the molecule using coupled-decoupled configurational-bias without requiring the backbone of the molecule to be aligned with the cavity created by the molecule that was removed.

The efficiency of the various molecular exchange algorithms is shown in Fig. 9 as a function of Monte Carlo step for x_butane = 0.1, 0.5, and 0.9. Uncertainties shown are the average over uncertainties for each histogram bin in the probability distribution. Both the ME-2 and ME-3 algorithms show that convergence of the probability distributions was achieved within 10 × 10⁶ MCS, while for ME-1 and configurational-bias insertions, convergence was not achieved within 20 × 10⁶ MCS. Depending on the composition, ME-3 provides efficiencies that are between 12 and 200 times greater than configurational-bias insertions for the insertion of perfluorobutane. Based on the trajectory of the uncertainties, it is unlikely that convergence of the probability distributions using standard Monte Carlo insertions would ever occur. Despite the fact that the ME-2 method provides slightly better acceptance rates for the molecular exchange move, at most compositions, ME-3 produces slightly faster convergence and better efficiencies. By growing the inserted molecule with coupled-decoupled configurational-bias,¹⁴ larger rearrangements take place in the system, even though more of the trial moves are rejected than in ME-2.

In Fig. 10, the probability distributions resulting from GCMC simulations with the various ME methods are presented for x_butane = 0.5, while data for x_butane = 0.1 and 0.9 are given in Figs. S10 and S11 of the supplementary material. The figure shows rapid convergence of the probability distributions for the ME-2 and ME-3 methods, while ME-1 and standard GCMC have not converged in 20 × 10⁶ MCS, although the uncertainties calculated for ME-1 are approximately half those of standard GCMC. In Fig. 11, heat maps are presented for the particle numbers and potential energies sampled during a liquid phase GCMC simulation. The heat maps illustrate how simulations with only configurational-bias insertions/deletions may become trapped in metastable states, resulting in poor sampling. Inclusion of the ME-3 algorithm produced a short equilibration period and a much broader sampling of the N₁, N₂, E phase space.

In order to compare the performance of the MEMC move with other advanced sampling techniques, such as the CBMC swap+identity switch (S+IS),³⁶ continuous fractional component Monte Carlo (CFCMC),^42,43 and configurational-bias continuous fractional component Monte Carlo (CB-CFCMC),⁴² the vapor-liquid coexistence curve for SPC/E water⁶² was predicted from the critical temperature to 0.44T_c. To enhance the acceptance rate for insertions and deletions of water and to provide a uniform basis for comparison, the strategy of Bai and Siepmann was used.³⁶ For regular CBMC swaps, oxygen was inserted first, followed by the two hydrogen atoms. 16 trials were used for the first atom and 8 trials for all remaining atoms. Simulations were performed as a mixture that contained approximately 0-10 “impurity” molecules, where the impurity molecule had an identical geometry to the SPC/E water model, but with partial charges reduced by a factor of 2 and the oxygen atom Lennard-Jones epsilon reduced by a factor of 4 compared to SPC/E water. Swap moves were performed only for impurity molecules, while the MEMC move was performed to exchange the impurity with water and vice versa. Move frequencies were adjusted to yield approximately to the same number of accepted molecule transfers for the swap and MEMC moves. Due to the poor performance of the ME-1 method in prior calculations, only the performance of the ME-2 and ME-3 methods was evaluated. An exchange ratio of one to one was used for all calculations.

The phase diagram for SPC/E water predicted from GCMC simulations using the ME-2 or ME-3 algorithm is shown in Fig. 12, with a comparison to prior simulations.⁶³ Additional data for the vapor pressure is provided in Fig. S12 of the supplementary material. Excellent agreement was observed, validating both the MEMC algorithms and the simulation code used to perform the calculations.

To compare the performance of MEMC with other methods, the effective number of molecule transfers was calculated. The effective number of molecule transfers was defined as the insertion of an impurity molecule by the swap move and its conversion to a regular water molecule by the MEMC move, or the conversion of regular water to impurity via MEMC and then deletion of impurity by the swap move. Exchanges of impurity to water and back to impurity were not counted. The effective acceptance rate was calculated from the effective number of molecule transfers divided by the sum of attempted swap and MEMC moves. The results of these calculations are summarized in Table II, with comparisons to the work of Bai and Siepmann³⁶ and Torres-Knoop et al.⁴² At 283 K, the effective acceptance rates for the ME-2 and ME-3 algorithms are 7.6 and 1.4 times greater, respectively, than the S+IS algorithm.³⁶ While the S+IS method reuses atomic coordinates of the molecule to be removed, the MEMC methods perform multiple trial orientations to insert the water molecule. In ME-2, first, the center of sub-volume was placed at the center of mass of the impurity, second, the z-axis of sub-volume was aligned with O–H bond of impurity, and then multiple rotational trials were performed around the z-axis of sub-volume. Aligning the O–H bond of water and sub-volume allows some of the original hydrogen bonding to be maintained, while finding an energetically favorable position for the oxygen atom through rotational trials around the z-axis of sub-volume, leading to significant improvements in the effective acceptance. In the ME-3 method, the oxygen atom of water was placed at the center of mass of the impurity molecule, and multiple rotational trials were performed on the sphere to find the most energetically favorable position. In order to maintain the hydrogen bonding formed by the impurity molecule, a large number of rotational trials are required, leading to a significant decrease in the acceptance efficiency compared to the ME-2 method.

Compared to the original CFCMC method of Shi and Maginn,⁴³ at 280 K, the ME-2 method exhibits twice the effective acceptance rate, while the ME-3 method is approximately 40% lower. The continuous fractional component Monte Carlo (CFCMC) and configurational-bias continuous fractional component Monte Carlo (CB-CFCMC) methods of Torres-Knoop et al.⁴² produced the largest acceptance rates of all methods. At 280 K, CFCMC and CB-CFCMC had acceptance rates that were 2.25 and 3.6 times larger, respectively, than the ME-2 method.

The acceptance efficiency was defined as the effective number of molecules transferred, divided by the total CPU time spent on swap and MEMC moves. In order to have a fair comparison between the acceptance efficiency of MEMC and S+IS, CFCMC, and CB-CFCMC methods, this quantity was normalized with respect to the acceptance efficiency of the standard CBMC method, minimizing the impact of the CPU choice on the relative performance of the algorithms. The results of these calculations are listed in Table III. At 280 K, the ME-2 method outperformed S+IS by 3.8 times, while the S+IS method is 23.9% better than ME-3. The performance of CFCMC and CB-CFCMC is 5-6 times greater than ME-2, although it should be noted that the acceptance rates reported by Torres-Knoop et al. for standard swaps of water were approximately 2.4 times lower than those reported in this work or in the work of Bai and Siepmann.³⁶ 

As mentioned earlier, achieving a statistically valid number of molecule insertions in low temperature $T<0.7Tc$ simulations of branched alkanes can be challenging. Here, 2,2,4-trimethylpentane is used as an example to highlight how the MEMC move can significantly extend the range of temperatures where GCMC simulations may be used to predict the vapor-liquid coexistence for a highly branched molecule. In this case, neopentane is used as the impurity molecule based on its similar structure to part of 2,2,4-trimethylpentane. This also illustrates the general nature of the MEMC algorithm, which does not require the molecules to be exchanged to be integer numbers of each other. In Fig. 13, the vapor-liquid coexistence curve for 2,2,4-trimethylpentane, using the ME-2 algorithm and GCMC+histogram reweighting Monte Carlo simulations, is presented. Additional data for the ME-3 algorithms are presented in Fig. S13 of the supplementary material. Using the ME-2 or ME-3 algorithms, it is possible to predict vapor-liquid coexistence to 280 K (0.51T_c), while prior simulations using only coupled-decoupled configurational bias were limited to 390 K (0.7T_c). In Table IV, a detailed comparison is presented for the acceptance rates for direct swaps of neopentane and 2,2,4-trimethylpentane, MEMC moves, effective acceptance rates, and effective acceptance rates per CPU time. Effective acceptance rate and acceptance efficiency were calculated using a similar method explained in Sec. IV C. The results of additional calculations performed with different CBMC parameters are given in Table S6 of the supplementary material. At all temperatures, the combination of the impurity swap plus the ME-2 or ME-3 method outperforms standard configurational-bias Monte Carlo. At 280 K, the relative acceptance efficiency (impurity swap+MEMC/standard CBMC) was 409 for ME-2 and 154 for ME-3. The ME-2 is more effective than ME-3 for branched molecules because it inserts the entire molecule at the same time and aligns the backbone of the molecule to be inserted with the backbone of the molecule to be removed. ME-3 regrows the entire molecule using coupled-decoupled CBMC; however, many of these growths fail because they are unable to satisfy the internal molecular constraints due to the bond bending and torsional potentials.⁶⁵ In future work, it may be possible to improve the performance of the ME-3 algorithm for branched molecules by inclusion of the Jacobian-Gaussian scheme⁶⁶ for generating bending angle trials in the CBMC growth.

In this work, three variants of the molecular exchange method were developed, which could be used to evaluate the efficiency of various aspects of the algorithms. Locating the exchange sub-volume randomly (ME-1) was found to have the lowest efficiency, since frequently no small molecules were found in the sub-volume that could be used for the molecular exchange, resulting in an immediate rejection of the move. The ME-1 method is suitable only for systems that are very dilute with respect to the concentration of the large molecule. By identifying a small molecule at random, first, placing the center of the sub-volume at the center of mass of the small molecule (ME-2) and aligning the backbone of the large molecule to be inserted with the small molecule to be removed, acceptance rates for the exchange move increased substantially. For water, the acceptance efficiency of the ME-2 method was found to be nearly 40 times greater than standard configurational-bias insertions, while for 2,2,4-trimethylpentane a 409 times improvement in acceptance efficiency was achieved. In the latter case, this was due to the use of a rigid-body insertion in ME-2, which eliminated the need to regrow the molecule in place. Finally, the inclusion of coupled-decoupled configurational-bias methods¹⁴ to grow sections of the molecule from a central atom (ME-3) placed at the center of the sub-volume resulted in a significant improvement in statistical efficiency compared to standard configurational-bias insertions for linear molecules without strong directional interactions. Improvements in efficiency of up to 200 times were observed for the perfluorobutane+n-butane system.

The algorithms presented in this work are notable because they were designed to work for any molecular topology over a wide range of compositions. Substantial performance gains were observed for ME-2 and ME-3 for all systems and compositions studied. As shown through the various case studies, however, each method has its strengths and weaknesses. For linear, non-polar molecules, ME-3 is generally more efficient than ME-2, while ME-2 offers better performance for small polar molecules, such as water, and highly branched molecules. Each of the algorithms was implemented, and the algorithms are now available in the open-source Monte Carlo simulation engine GOMC, which is available to the public at GitHub.³⁷ 

See supplementary material for additional details regarding the implementation of the MEMC moves, force field parameters, and tabulated numerical data.

We would like to thank Ariana Torres-Knoop and Thijs Vlugt for sharing tabulated numerical data from their work on optimizing the performance of continuous fractional component Monte Carlo and Peng Bai for helpful discussions. Funding from the National Science Foundation No. OAC-1642406 is gratefully acknowledged. Some of the computations in this work were performed with resources from the Grid Computing initiative at Wayne State University.

The authors declare no competing financial interest.