A generalized identity exchange algorithm is presented for Monte Carlo simulations in the grand canonical ensemble. The algorithm, referred to as molecular exchange Monte Carlo, may be applied to multicomponent systems of arbitrary molecular topology and provides significant enhancements in the sampling of phase space over a wide range of compositions and temperatures. Three different approaches are presented for the insertion of large molecules, and the pros and cons of each method are discussed. The performance of the algorithms is highlighted through grand canonical Monte Carlo histogram-reweighting simulations performed on a number of systems, which include methane+n-alkanes, butane+perfluorobutane, water+impurity, and 2,2,4-trimethylpentane+neopentane. Relative acceptance efficiencies for molecule transfers of up to 400 times that of standard configurational-bias Monte Carlo are obtained.

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 mixtures9,10 and have been used to determine critical micelle concentrations for model surfactants.11 

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.12 The introduction of configurational-bias Monte Carlo (CBMC) enabled the successful simulation of chain molecule adsorption in zeolites,13 which was followed by the coupled-decoupled14 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.7Tc were approximately 0.3%.9 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.24 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.31 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,26 a polymer composed of solvent monomers,31 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 NL large and NS small molecules, an “insertion move” is an attempt to exchange one large molecule with NEX small molecules inside a predefined exchange sub-volume VEX, and a “deletion move” is an attempt to exchange NEX 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 NEX and VEX “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

Koldnew=Knewold,
(1)

where K(ij) 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

Koldnew=Nold×αoldnew×accoldnew.
(2)

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

accoldnewaccnewold=NnewNold×αnewoldαoldnew.
(3)

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

NnewNold=eβUneweβμLNL1+μSNS+NEXeβUoldeβNLμL+NSμS=eβNEXμSμLeβUnewUold,
(4)

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 NEX 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

NnewNold=eβUneweβμLNL+1+μSNSNEXeβUoldeβNLμL+NSμS=eβμLNEXμSeβUnewUold.
(5)

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 NEX particular small molecules, the probability of choosing a particular large molecule, the probability of generating trial configurations for NEX small molecules, and the probability of generating trial configurations for the large molecule,

αoldnew=Psubvoldnew×PpickSoldnew×PpickLoldnew×PposSoldnew×PposLoldnew.
(6)

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 VEX 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 VEX 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.

FIG. 1.

Schematic of the ME-1 algorithm. Selected or inserted molecule (green), trial position (light red), and actual position of the molecule (solid red). Top row represents the exchange of two small molecules with one large molecule (insertion). The exchange sub-volume is defined as the orange box. (a) Identifying small molecules within the sub-volume with a random geometric center and orientation. (b) Generating CBMC trials (rotation and COM location) for one of the small molecules and then removing it. (c) Generating CBMC trials (rotation and COM location) for the second small molecule and then removing it. (d) Aligning the backbone of the large molecule with the sub-volume and performing CBMC rotational trials around the z-axis of the sub-volume. Bottom row represents the exchange of a large molecule (deletion) with two small molecules. (a) Aligning the sub-volume with the large molecule’s backbone with the geometric center placed at COM of the large molecule and identifying the small molecules within the sub-volume. (b) Generating CBMC rotational trials for the large molecule around the z-axis of the sub-volume and then removing it. (c) Generating CBMC trials (rotation and COM location) for the first small molecule and then inserting it into the sub-volume. (d) Generating CBMC trials (rotation and COM location) for the second small molecule and then inserting it into the sub-volume.

FIG. 1.

Schematic of the ME-1 algorithm. Selected or inserted molecule (green), trial position (light red), and actual position of the molecule (solid red). Top row represents the exchange of two small molecules with one large molecule (insertion). The exchange sub-volume is defined as the orange box. (a) Identifying small molecules within the sub-volume with a random geometric center and orientation. (b) Generating CBMC trials (rotation and COM location) for one of the small molecules and then removing it. (c) Generating CBMC trials (rotation and COM location) for the second small molecule and then removing it. (d) Aligning the backbone of the large molecule with the sub-volume and performing CBMC rotational trials around the z-axis of the sub-volume. Bottom row represents the exchange of a large molecule (deletion) with two small molecules. (a) Aligning the sub-volume with the large molecule’s backbone with the geometric center placed at COM of the large molecule and identifying the small molecules within the sub-volume. (b) Generating CBMC rotational trials for the large molecule around the z-axis of the sub-volume and then removing it. (c) Generating CBMC trials (rotation and COM location) for the first small molecule and then inserting it into the sub-volume. (d) Generating CBMC trials (rotation and COM location) for the second small molecule and then inserting it into the sub-volume.

Close modal

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 VEX, with its geometric center located randomly within the simulation box of volume V (with the probability proportional to V−1) and a random orientation. Determine the total number of small molecules within the exchange sub-volume (NS,VEX) based on their center of mass.

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

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

  4. Repeat the steps a and b for NEX cycles (i = 1, 2,…, NEX) 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 VEX. 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=1jr=1kexpβUi,p,r, where Ui,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 Ui,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. Pi,old is the probability of inserting the ith small molecule back in its original configuration in the reverse move (newold).

  5. Insert the COM of the large molecule at the geometric center of the exchange sub-volume VEX 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=1kexpβUr, where Ur 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 NL large molecules within the simulation box with the probability of 1/NL.

  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 NS,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=1kexpβUr, where Ur 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 Uk is the interaction energy of the large molecule at the original orientation with all other molecules in the simulation box. Pold is the probability of inserting the large molecule at its original configuration in the reverse move (newold). Then remove the large molecule from the simulation box.

  4. Repeat the steps a → c for NEX cycles (i = 1, 2,…, NEX) 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 VEX.

    • 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=1jr=1kexpβUi,p,r, where Ui,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 (NL + 1, NSNEX) → old (NL, NS) to that of the reverse move is

αnewoldαoldnew=1NL+11V×NEX!VEXNEXNEX!NS,VEXNEX!NS,VEX!×i=1NEXPi,oldPnew.
(7)

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

accoldnew=min1,VNL+1×NS,VEX!VEXNEXNS,VEXNEX!×Wnewi=1NEXWi,old×eβμLNEXμS.
(8)

For the large molecule deletion move, the ratio of the probability of generating the move new (NL − 1, NS + NEX) → old (NL, NS) to that of the reverse move is

αnewoldαoldnew=1V1NL×NEX!NS,VEX!NS,VEX+NEX!NEX!VEXNEX×i=1NEXPi,newPold.
(9)

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

accoldnew=min1,NLV×VEXNEX×NS,VEX!NS,VEX+NEX!×i=1NEXWi,newWold×eβNEXμSμL.
(10)

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.30 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 NEX 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 VEX 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 VEX is placed on the COM of a randomly selected small molecule. If the small molecule is monoatomic, the orientation of VEX 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.

FIG. 2.

Schematic of the ME-2 algorithm. Selected or inserted molecule (green), trial position (light red), and actual position of the molecule (solid red). Top row represents the exchange of two small molecules with one large molecule (insertion). The sub-volume is defined as the orange box. (a) Aligning the sub-volume with a randomly selected small molecule’s backbone with the geometric center placed at COM of the selected small molecule and identifying the small molecules within the sub-volume. (b) Generating CBMC trials (rotation and COM location) for one of the small molecules and then removing it. (c) Generating CBMC rotational trials around the selected small molecule and then removing it. (d) Aligning the backbone of the large molecule with the sub-volume and performing CBMC rotational trials around the z-axis of the sub-volume. Bottom row represents the exchange of one large molecule with two small molecules (deletion). (a) Aligning the sub-volume with the large molecule’s backbone with the geometric center placed at COM of the large molecule and identifying the small molecules within the sub-volume. (b) Generating CBMC rotational trials for the large molecule around the z-axis of the sub-volume and then removing it. (c) Placing the COM of the first small molecule at the geometric center of the sub-volume and generating the CBMC rotational trials around the z-axis and then inserting it into the sub-volume. (d) Generating CBMC trials (rotation and COM location) for the second small molecule and then inserting it into the sub-volume.

FIG. 2.

Schematic of the ME-2 algorithm. Selected or inserted molecule (green), trial position (light red), and actual position of the molecule (solid red). Top row represents the exchange of two small molecules with one large molecule (insertion). The sub-volume is defined as the orange box. (a) Aligning the sub-volume with a randomly selected small molecule’s backbone with the geometric center placed at COM of the selected small molecule and identifying the small molecules within the sub-volume. (b) Generating CBMC trials (rotation and COM location) for one of the small molecules and then removing it. (c) Generating CBMC rotational trials around the selected small molecule and then removing it. (d) Aligning the backbone of the large molecule with the sub-volume and performing CBMC rotational trials around the z-axis of the sub-volume. Bottom row represents the exchange of one large molecule with two small molecules (deletion). (a) Aligning the sub-volume with the large molecule’s backbone with the geometric center placed at COM of the large molecule and identifying the small molecules within the sub-volume. (b) Generating CBMC rotational trials for the large molecule around the z-axis of the sub-volume and then removing it. (c) Placing the COM of the first small molecule at the geometric center of the sub-volume and generating the CBMC rotational trials around the z-axis and then inserting it into the sub-volume. (d) Generating CBMC trials (rotation and COM location) for the second small molecule and then inserting it into the sub-volume.

Close modal

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 NS small molecules in the simulation box with the probability of 1/NS. This molecule will be the last molecule to be removed from the system.

  2. Define VEX 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 VEX is assigned randomly. Determine the number of small molecules NS,VEX within VEX (NS,VEX includes the molecules selected in step 1).

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

  4. Select NEX − 1 small molecules out of NS,VEX − 1, with the probability (NEX1)!(NS,VEXNEX)!/(NS,VEX1)!.

  5. Repeat steps a and b of the large molecule insertion move of ME-1 for NEX − 1 cycles (i = 1, 2,…, NEX − 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=1kexpβUNEX,r, where UNEX,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 UNEX,k is the interaction energy of the last small molecule at its original configuration with all other molecules remaining in the simulation box. Pi,old is the probability of inserting the ith small molecule back in its original configuration in the reverse move (newold).

  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 VEX 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=1kexpβU1,r, where U1,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 NEX − 1 cycles (i = 2,…, NEX) to insert the small molecules with the probability NEX1!/VEXNEX1.

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

αnewoldαoldnew=1NL+11NS×NEX1!VEXNEX1NEX1!NS,VEXNEX!NS,VEX1!×i=1NEXPi,oldPnew.
(11)

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

acc(oldnew)=min1,NSNL+1×(NS,VEX1)!VEX(NEX1)(NS,VEXNEX)!×Wnewi=1NEXWi,old×eβ[μLNEXμS].
(12)

For the large molecule deletion move, the ratio of the probability of generating configuration newNL1,NS+NEXoldNL,NS to that of the reverse move is

αnewoldαoldnew=1NS+NEX1NL×NEX1!NS,VEX!NS,VEX+NEX1!NEX1!VEXNEX1×i=1NEXPi,newPold.
(13)

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

acc(oldnew)=min1,NL(NS+NEX)×VEX(NEX1)×NS,VEX!(NS,VEX+NEX1)!×i=1NEXWi,newWold×eβ[NEXμSμL].
(14)

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

acc(NLNL+1)=min1,NS(NL+1)×WnewWold×eβ[μLμS],
(15)
acc(NLNL1)=min1,NL(NS+1)×WnewWold×eβ[μSμL].
(16)

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 VEX. 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 VEX and the molecule is built segment by segment using the coupled-decoupled configurational-bias Monte Carlo (CBMC) algorithm.14 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 Wold of the large molecule is calculated. Insertion and deletion of NEX 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 NS small molecules in the simulation box with the probability 1/NS. This molecule will be the last molecule to be removed from the system.

  2. Define an orthogonal exchange sub-volume VEX 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 NS,VEX within VEX (NS,VEX includes the molecule selected in step 1).

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

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

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

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 NL with the probability of 1/NL.

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

  3. Perform coupled-decoupled CBMC for the large molecule and calculate the Rosenbluth weight Wold and Pold.

  4. Repeat steps 2-4 of ME-2 to insert NEX small molecules within VEX.

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).

FIG. 3.

Schematic of the ME-3. Selected or inserted molecule (green), trial position (light red), and actual position of the molecule (solid red). Top row represents the exchange of two small molecules with one large molecule (insertion). The sub-volume is defined as the orange box. (a) Defining the sub-volume with a random orientation, where its geometric center is placed at a randomly selected small molecule’s COM and identifying the small molecules within the sub-volume. (b) Generating CBMC trials (rotation and COM location) for one of the small molecules and then removing it. (c) Generating CBMC rotational trials around its COM of the selected small molecule and then removing it. (d) Placing the predefined atom of the large molecule at the geometric center of the sub-volume and growing the large molecule using the coupled-decoupled CBMC technique. Bottom row represents the exchange of a large molecule with two small molecules (deletion). (a) Defining the sub-volume with a random orientation with the geometric center placed at the predefined atom of the large molecule and identifying the small molecules within the sub-volume. (b) Generating coupled-decoupled CBMC trials and then removing it. (c) Placing the COM of the first small molecule at the geometric center of the sub-volume, generating CBMC rotational trials around its COM, and then inserting it into the sub-volume. (d) Generating CBMC trials (rotation and COM location) for the second small molecule and then inserting it into the sub-volume.

FIG. 3.

Schematic of the ME-3. Selected or inserted molecule (green), trial position (light red), and actual position of the molecule (solid red). Top row represents the exchange of two small molecules with one large molecule (insertion). The sub-volume is defined as the orange box. (a) Defining the sub-volume with a random orientation, where its geometric center is placed at a randomly selected small molecule’s COM and identifying the small molecules within the sub-volume. (b) Generating CBMC trials (rotation and COM location) for one of the small molecules and then removing it. (c) Generating CBMC rotational trials around its COM of the selected small molecule and then removing it. (d) Placing the predefined atom of the large molecule at the geometric center of the sub-volume and growing the large molecule using the coupled-decoupled CBMC technique. Bottom row represents the exchange of a large molecule with two small molecules (deletion). (a) Defining the sub-volume with a random orientation with the geometric center placed at the predefined atom of the large molecule and identifying the small molecules within the sub-volume. (b) Generating coupled-decoupled CBMC trials and then removing it. (c) Placing the COM of the first small molecule at the geometric center of the sub-volume, generating CBMC rotational trials around its COM, and then inserting it into the sub-volume. (d) Generating CBMC trials (rotation and COM location) for the second small molecule and then inserting it into the sub-volume.

Close modal

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.37 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.38 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.37 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.39 Psfgen was used to generate coordinate (*.pdb) and connectivity (*.psf) files.40 Potentials were truncated at 10 Å, and analytical tail corrections were applied.41 To enhance the acceptance rate for molecule insertions, the coupled-decoupled configurational-bias Monte Carlo (CBMC) algorithm was used.14 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 × 107 Monte Carlo steps (MCS), without equilibration period. Simulations used to generate phase diagrams were run for 5 × 107 MCS with a 5 × 106 MCS equilibration period. Every 200-500 MCS, the instantaneous state of the system (N1, N2, E) was saved as a histogram. Every one million MCS, the natural log of distribution of large particles ln(PN) 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,

η=σ2s1,
(17)

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.7Tc. For pure fluid phase diagrams, calculations were performed from the critical temperature to 0.44Tc – 0.51Tc. Performing grand canonical Monte Carlo simulations, using standard configurational-bias methods,14 below 0.7Tc 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.7Tc can be done using standard configurational-bias methods in grand canonical or Gibbs ensemble Monte Carlo simulations.9,52–54 However, below 0.7Tc, 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.

FIG. 4.

Probability distributions predicted from gas (μbutane = −2960, μmethane = −2000) and liquid (μbutane = −2840, μmethane = −2000) phase GCMC simulations of methane+n-butane at 277 K. Solid lines denote the probability distributions for n-butane (black) and methane (blue) using standard configurational-bias insertions and deletions. Dashed lines denote the probability distributions for n-butane (red) and methane (green) using the ME-3 algorithm.

FIG. 4.

Probability distributions predicted from gas (μbutane = −2960, μmethane = −2000) and liquid (μbutane = −2840, μmethane = −2000) phase GCMC simulations of methane+n-butane at 277 K. Solid lines denote the probability distributions for n-butane (black) and methane (blue) using standard configurational-bias insertions and deletions. Dashed lines denote the probability distributions for n-butane (red) and methane (green) using the ME-3 algorithm.

Close modal

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

FIG. 5.

Pressure composition diagram for methane+n-butane at 277 K predicted from GCMC+histogram reweighting simulations using Mie potentials.52 Experimental data (circles),51 standard configurational-bias insertions (red lines), ME-3 algorithm (green lines).

FIG. 5.

Pressure composition diagram for methane+n-butane at 277 K predicted from GCMC+histogram reweighting simulations using Mie potentials.52 Experimental data (circles),51 standard configurational-bias insertions (red lines), ME-3 algorithm (green lines).

Close modal

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 xmethane = 0.7 to 70× for methane+ethane at xmethane = 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.

TABLE I.

n-alkane insertion/removal acceptance percentages in GCMC liquid phase simulations of methane+n-alkane mixtures for CBMC, ME-1, ME-2, and ME-3 methods.

Binary systemSub-volume size (Å)NEXxCH4CBMCME-1ME-2ME-3
methane+ethane 5 × 5 × 6 0.1 0.33 0.11 11.68 23.62 
   0.5 1.47 0.96 16.20 33.33 
   0.9 8.3 4.18 24.09 47.84 
methane+propane 5 × 5 × 7 0.1 0.08 0.05 3.42 4.13 
   0.4 0.38 0.40 5.67 7.21 
   0.8 5.18 3.36 13.56 18.36 
methane+n-butane 5 × 5 × 8 0.1 0.14 0.025 0.835 2.373 
   0.3 0.33 0.099 1.207 3.421 
   0.6 2.52 0.948 3.378 8.128 
 5 × 5 × 8 0.1 0.14 0.019 0.196 0.362 
   0.3 0.33 0.144 0.557 0.928 
   0.6 2.52 1.262 2.288 3.160 
 8.8 × 8.8 × 11.8 0.1 0.14 0.022 0.398 0.984 
   0.3 0.33 0.086 0.821 1.860 
   0.6 2.52 0.621 2.682 5.252 
methane+n-pentane 5 × 5 × 9 0.1 0.064 0.007 0.209 0.824 
   0.5 0.397 0.116 0.638 2.163 
   0.7 2.461 0.666 1.72 4.814 
 5 × 5 × 9 0.1 0.064 0.006 0.086 0.189 
   0.5 0.397 0.270 0.736 1.160 
   0.7 2.461 1.332 2.389 3.170 
 8.8 × 8.8 × 13 0.1 0.064 0.008 0.145 0.455 
   0.5 0.397 0.102 0.675 1.806 
   0.7 2.461 0.473 2.054 4.133 
Binary systemSub-volume size (Å)NEXxCH4CBMCME-1ME-2ME-3
methane+ethane 5 × 5 × 6 0.1 0.33 0.11 11.68 23.62 
   0.5 1.47 0.96 16.20 33.33 
   0.9 8.3 4.18 24.09 47.84 
methane+propane 5 × 5 × 7 0.1 0.08 0.05 3.42 4.13 
   0.4 0.38 0.40 5.67 7.21 
   0.8 5.18 3.36 13.56 18.36 
methane+n-butane 5 × 5 × 8 0.1 0.14 0.025 0.835 2.373 
   0.3 0.33 0.099 1.207 3.421 
   0.6 2.52 0.948 3.378 8.128 
 5 × 5 × 8 0.1 0.14 0.019 0.196 0.362 
   0.3 0.33 0.144 0.557 0.928 
   0.6 2.52 1.262 2.288 3.160 
 8.8 × 8.8 × 11.8 0.1 0.14 0.022 0.398 0.984 
   0.3 0.33 0.086 0.821 1.860 
   0.6 2.52 0.621 2.682 5.252 
methane+n-pentane 5 × 5 × 9 0.1 0.064 0.007 0.209 0.824 
   0.5 0.397 0.116 0.638 2.163 
   0.7 2.461 0.666 1.72 4.814 
 5 × 5 × 9 0.1 0.064 0.006 0.086 0.189 
   0.5 0.397 0.270 0.736 1.160 
   0.7 2.461 1.332 2.389 3.170 
 8.8 × 8.8 × 13 0.1 0.064 0.008 0.145 0.455 
   0.5 0.397 0.102 0.675 1.806 
   0.7 2.461 0.473 2.054 4.133 

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 (NEX = 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 × 107 MCS, with average uncertainties of approximately 0.05. The ME-1 algorithm and configurational-bias insertions show similar convergence properties. However, with 2 × 107 MCS, each produced uncertainties that were approximately double than those of the ME-3 and ME-2 methods.

FIG. 6.

Efficiency and standard deviation in methane+n-butane at 255 K. Lines represent the efficiency and average uncertainty in probability distributions generated from GCMC simulations. Standard configurational-bias insertions (black), ME-1 (red), ME-2 (green), and ME-3 (blue). The MEMC move was performed with the exchange ratio of one butane with one methane.

FIG. 6.

Efficiency and standard deviation in methane+n-butane at 255 K. Lines represent the efficiency and average uncertainty in probability distributions generated from GCMC simulations. Standard configurational-bias insertions (black), ME-1 (red), ME-2 (green), and ME-3 (blue). The MEMC move was performed with the exchange ratio of one butane with one methane.

Close modal

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 xmethane = 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 × 106 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.

FIG. 7.

Probability distributions for methane+n-butane at 255 K and xmethane = 0.3. After simulations of 1 × 106 MCS (magenta), 5 × 106 MCS (green), 1 × 107 MCS (blue), 1.5 × 107 MCS (red), and 2 × 107 MCS (black), (a) standard configurational-bias insertions, (b) ME-1, (c) ME-2, and (d) ME-3.

FIG. 7.

Probability distributions for methane+n-butane at 255 K and xmethane = 0.3. After simulations of 1 × 106 MCS (magenta), 5 × 106 MCS (green), 1 × 107 MCS (blue), 1.5 × 107 MCS (red), and 2 × 107 MCS (black), (a) standard configurational-bias insertions, (b) ME-1, (c) ME-2, and (d) ME-3.

Close modal

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 C4F10 and 272.61 K for C4H10) 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),55 Perturbed Chain Statistical Associating Fluid Theory (PC-SAFT),56 and Group Contribution Statistical Associating Fluid Theory Variable Range (GC-SAFT-VR),57 which showed close agreement with experimental data.58 Gibbs ensemble Monte Carlo simulations using an identity exchange move have been used to study liquid-liquid equilibria for n-heptane+perfluoheptane,59 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,52 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 rules60,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.

FIG. 8.

Pressure-composition diagram for perfluorobutane+n-butane at 259.95 K. The predictions from GCMC+histogram reweighting simulations using the ME-2 algorithm are given by red lines, while experimental data58 are represented by black circles. The line connecting the experimental data points is provided as guide to the eye.

FIG. 8.

Pressure-composition diagram for perfluorobutane+n-butane at 259.95 K. The predictions from GCMC+histogram reweighting simulations using the ME-2 algorithm are given by red lines, while experimental data58 are represented by black circles. The line connecting the experimental data points is provided as guide to the eye.

Close modal

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,14 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 xbutane = 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 × 106 MCS, while for ME-1 and configurational-bias insertions, convergence was not achieved within 20 × 106 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,14 larger rearrangements take place in the system, even though more of the trial moves are rejected than in ME-2.

FIG. 9.

Efficiency and standard deviation in the perfluorobutane+n-butane binary mixture at 259.95 K. Lines represent the efficiency and average uncertainty in the perfluorobutane probability distribution; standard configurational-bias insertions (black), ME-1 (red), ME-2 (green), and ME-3 (blue). The MEMC moves were performed with an exchange ratio of one to one.

FIG. 9.

Efficiency and standard deviation in the perfluorobutane+n-butane binary mixture at 259.95 K. Lines represent the efficiency and average uncertainty in the perfluorobutane probability distribution; standard configurational-bias insertions (black), ME-1 (red), ME-2 (green), and ME-3 (blue). The MEMC moves were performed with an exchange ratio of one to one.

Close modal

In Fig. 10, the probability distributions resulting from GCMC simulations with the various ME methods are presented for xbutane = 0.5, while data for xbutane = 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 × 106 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 N1, N2, E phase space.

FIG. 10.

Molecule probability distribution for perfluorobutane+n-butane at xbutane = 0.5 and 259.95 K. After simulations of 1 × 106 MCS (magenta), 5 × 106 MCS (green), 1 × 107 MCS (blue), 1.5 × 107 MCS (red), and 2 × 107 MCS (black), (a) standard configurational-bias insertions, (b) ME-1, (c) ME-2, and (d) ME-3.

FIG. 10.

Molecule probability distribution for perfluorobutane+n-butane at xbutane = 0.5 and 259.95 K. After simulations of 1 × 106 MCS (magenta), 5 × 106 MCS (green), 1 × 107 MCS (blue), 1.5 × 107 MCS (red), and 2 × 107 MCS (black), (a) standard configurational-bias insertions, (b) ME-1, (c) ME-2, and (d) ME-3.

Close modal
FIG. 11.

Heat maps of particle numbers (left panel) and potential energies (right panel) sampled during liquid phase grand canonical Monte Carlo simulations of perfluorobutane+n-butane at 259.95 K. Upper figures correspond to GCMC simulations with standard configurational-bias insertions/deletions, while the bottom figures correspond to GCMC simulations with the ME-3 algorithm.

FIG. 11.

Heat maps of particle numbers (left panel) and potential energies (right panel) sampled during liquid phase grand canonical Monte Carlo simulations of perfluorobutane+n-butane at 259.95 K. Upper figures correspond to GCMC simulations with standard configurational-bias insertions/deletions, while the bottom figures correspond to GCMC simulations with the ME-3 algorithm.

Close modal

In order to compare the performance of the MEMC move with other advanced sampling techniques, such as the CBMC swap+identity switch (S+IS),36 continuous fractional component Monte Carlo (CFCMC),42,43 and configurational-bias continuous fractional component Monte Carlo (CB-CFCMC),42 the vapor-liquid coexistence curve for SPC/E water62 was predicted from the critical temperature to 0.44Tc. 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.36 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.63 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.

FIG. 12.

Vapor-liquid coexistence curve for SPC/E water predicted from GCMC+histogram reweighting simulations. NIST Chemistry WebBook64 (solid lines), values obtained by Boulougouris et al.63 (green circles), ME-2 algorithm (red squares), and ME-3 algorithm (blue triangles).

FIG. 12.

Vapor-liquid coexistence curve for SPC/E water predicted from GCMC+histogram reweighting simulations. NIST Chemistry WebBook64 (solid lines), values obtained by Boulougouris et al.63 (green circles), ME-2 algorithm (red squares), and ME-3 algorithm (blue triangles).

Close modal

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 Siepmann36 and Torres-Knoop et al.42 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.36 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.

TABLE II.

Comparison of Swap+MEMC move acceptance percentages with standard CBMC, S+IS,36 CFCMC,42,43 and CB-CFCMC42 for SPC/E water.

%PImpacc%Pwateracc%Pwateracc%Pwateracc
(CBMC)%PSwitchacc%PEffectiveacc(CBMC)(CFCMC)(CBCFCMC)
BaiBaiTorres-KnoopShiTorres-KnoopTorres-Knoop
T (K)This worket al.ME-2ME-3ISME-2ME-3S+ISThis worket al.et al.et al.et al.et al.
280 5.7 … 5.70 0.59 … 2.73 0.51 … 0.063 … 0.027 1.38 6.16 9.86 
283 5.9 4.3 6.07 0.61 1.4 2.94 0.53 0.36 0.076 0.061 … … … … 
313 6.3 … 6.74 0.98 … 3.35 0.83 … 0.167 … 0.068 1.00 7.49 11.7 
343 6.8 7.8 6.61 1.10 3.1 3.28 0.91 0.73 0.35 0.37 … … … … 
348 7.0 … 6.47 1.28 … 2.94 1.07 … 0.423 … 0.155 2.18 9.52 14.93 
375 9.8 … 8.67 2.11 … 4.55 1.71 … 0.761 … 0.286 … 10.14 16.53 
473 20.5 22 14.84 6.31 7.3 8.48 4.84 2.2 3.989 3.5 1.374 1.98 15.17 21.82 
500 23 … 15.95 7.49 … 9.29 5.62 … 5.556 … 1.964 … 15.23 21.5 
%PImpacc%Pwateracc%Pwateracc%Pwateracc
(CBMC)%PSwitchacc%PEffectiveacc(CBMC)(CFCMC)(CBCFCMC)
BaiBaiTorres-KnoopShiTorres-KnoopTorres-Knoop
T (K)This worket al.ME-2ME-3ISME-2ME-3S+ISThis worket al.et al.et al.et al.et al.
280 5.7 … 5.70 0.59 … 2.73 0.51 … 0.063 … 0.027 1.38 6.16 9.86 
283 5.9 4.3 6.07 0.61 1.4 2.94 0.53 0.36 0.076 0.061 … … … … 
313 6.3 … 6.74 0.98 … 3.35 0.83 … 0.167 … 0.068 1.00 7.49 11.7 
343 6.8 7.8 6.61 1.10 3.1 3.28 0.91 0.73 0.35 0.37 … … … … 
348 7.0 … 6.47 1.28 … 2.94 1.07 … 0.423 … 0.155 2.18 9.52 14.93 
375 9.8 … 8.67 2.11 … 4.55 1.71 … 0.761 … 0.286 … 10.14 16.53 
473 20.5 22 14.84 6.31 7.3 8.48 4.84 2.2 3.989 3.5 1.374 1.98 15.17 21.82 
500 23 … 15.95 7.49 … 9.29 5.62 … 5.556 … 1.964 … 15.23 21.5 

Compared to the original CFCMC method of Shi and Maginn,43 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.42 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.36 

TABLE III.

Comparison of relative acceptance efficiency for the MEMC, S+IS,36 CFCMC,42 and CB-CFCMC42 methods.

T (K)ME-2ME-3S+ISCFCMC42 CB-CFCMC42 
280 38.8 7.61 … 243.47 195.28 
283 34.1 6.49 10 … … 
313 19.33 4.91 … 97.07 85.27 
343 11.04 3.32 3.45 … … 
348 7.97 3.02 … 52.18 42.69 
375 6.39 2.47 … 33.16 27.59 
473 2.08 1.25 1.23 7.74 6.85 
500 1.65 1.04 … 6.52 5.18 
T (K)ME-2ME-3S+ISCFCMC42 CB-CFCMC42 
280 38.8 7.61 … 243.47 195.28 
283 34.1 6.49 10 … … 
313 19.33 4.91 … 97.07 85.27 
343 11.04 3.32 3.45 … … 
348 7.97 3.02 … 52.18 42.69 
375 6.39 2.47 … 33.16 27.59 
473 2.08 1.25 1.23 7.74 6.85 
500 1.65 1.04 … 6.52 5.18 

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.51Tc), while prior simulations using only coupled-decoupled configurational bias were limited to 390 K (0.7Tc). 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.65 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 scheme66 for generating bending angle trials in the CBMC growth.

FIG. 13.

Vapor-liquid coexistence curve for 2,2,4-trimethylpentane predicted from GCMC+histogram reweighting simulations using Mie potentials.9 Experimental data (solid lines),67 ME-2 algorithm (red circles), and prior calculations using only coupled-decoupled configurational-bias Monte Carlo (green circles).9 

FIG. 13.

Vapor-liquid coexistence curve for 2,2,4-trimethylpentane predicted from GCMC+histogram reweighting simulations using Mie potentials.9 Experimental data (solid lines),67 ME-2 algorithm (red circles), and prior calculations using only coupled-decoupled configurational-bias Monte Carlo (green circles).9 

Close modal
TABLE IV.

Comparison of acceptance rates for swaps of the impurity molecule (neopentane), identity exchange via the MEMC algorithm, and swaps performed with standard configurational-bias Monte Carlo for 2,2,4-trimethylpentane.

Effective acceptanceRelative acceptance
%PImpacc%PSwitchacc%PEffectiveacc%Paccper CPU time (s−1)efficiency
T (K)SwapME-2ME-3ME-2ME-3CBMCCBMCME-2ME-3ME-2ME-3
280 0.013 0.89 0.03 0.013 0.008 0.00008 0.0003 0.109 0.041 409.2 153.7 
330 0.10 2.21 0.15 0.096 0.057 0.0008 0.0026 0.917 0.288 356.9 112.0 
390 0.85 5.69 0.55 0.653 0.274 0.022 0.0769 5.727 1.135 74.5 14.8 
450 4.09 9.84 1.27 2.645 0.837 0.225 0.838 24.12 3.497 28.8 4.17 
510 13.50 21.07 2.89 6.613 1.894 1.026 4.120 55.74 7.210 13.5 1.75 
Effective acceptanceRelative acceptance
%PImpacc%PSwitchacc%PEffectiveacc%Paccper CPU time (s−1)efficiency
T (K)SwapME-2ME-3ME-2ME-3CBMCCBMCME-2ME-3ME-2ME-3
280 0.013 0.89 0.03 0.013 0.008 0.00008 0.0003 0.109 0.041 409.2 153.7 
330 0.10 2.21 0.15 0.096 0.057 0.0008 0.0026 0.917 0.288 356.9 112.0 
390 0.85 5.69 0.55 0.653 0.274 0.022 0.0769 5.727 1.135 74.5 14.8 
450 4.09 9.84 1.27 2.645 0.837 0.225 0.838 24.12 3.497 28.8 4.17 
510 13.50 21.07 2.89 6.613 1.894 1.026 4.120 55.74 7.210 13.5 1.75 

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 methods14 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.37 

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.

1.
J. J.
de Pablo
,
M.
Laso
,
J. I.
Siepmann
, and
U. W.
Suter
,
Mol. Phys.
80
,
55
(
1993
).
2.
J. I.
Siepmann
and
D.
Frenkel
,
Mol. Phys.
75
,
59
(
1992
).
3.
C. E.
Wilmer
and
R. Q.
Snurr
,
Chem. Eng. J.
171
,
775
(
2011
).
4.
J.
Kim
,
L. C.
Lin
,
R. L.
Martin
,
J. A.
Swisher
,
M.
Haranczyk
, and
B.
Smit
,
Langmuir
28
,
11914
(
2012
).
5.
D.
Banerjee
,
C. M.
Simon
,
A. M.
Plonka
,
R. K.
Motkuri
,
J.
Liu
,
X. Y.
Chen
,
B.
Smit
,
J. B.
Parise
,
M.
Haranczyk
, and
P. K.
Thallapally
,
Nat. Commun.
7
,
11831
(
2016
).
6.
C. E.
Wilmer
,
M.
Leaf
,
C. Y.
Lee
,
O. K.
Farha
,
B. G.
Hauser
,
J. T.
Hupp
, and
R. Q.
Snurr
,
Nat. Chem.
4
,
83
(
2012
).
7.
A. M.
Ferrenberg
and
R. H.
Swendsen
,
Phys. Rev. Lett.
63
,
1195
(
1989
).
8.
A. M.
Ferrenberg
and
R. H.
Swendsen
,
Phys. Rev. Lett.
61
,
2635
(
1988
).
9.
J. R.
Mick
,
M. S.
Barhaghi
,
B.
Jackman
,
L.
Schwiebert
, and
J. J.
Potoff
,
J. Chem. Eng. Data
62
,
1806
(
2017
).
10.
J. R.
Errington
and
A. Z.
Panagiotopoulos
,
J. Phys. Chem. B
103
,
6314
(
1999
).
11.
M. A.
Floriano
,
E.
Caponetti
, and
A. Z.
Panagiotopoulos
,
Langmuir
15
,
3143
(
1999
).
12.
R. Q.
Snurr
,
A. T.
Bell
, and
D. N.
Theodorou
,
J. Phys. Chem.
97
,
13742
(
1993
).
14.
M. G.
Martin
and
J. I.
Siepmann
,
J. Phys. Chem. B
103
,
4508
(
1999
).
15.
J. R.
Errington
and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
111
,
9731
(
1999
).
16.
M. D.
Macedonia
and
E. J.
Maginn
,
Mol. Phys.
96
,
1375
(
1999
).
17.
M. R.
Stapleton
and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
92
,
1285
(
1990
).
18.
V.
Ortiz
,
J. R.
Maury-Evertsz
, and
G. E.
Lopez
,
Chem. Phys. Lett.
368
,
452
(
2003
).
20.
W.
Shi
and
E. J.
Maginn
,
J. Chem. Theory Comput.
3
,
1451
(
2007
).
21.
A.
Torres-Knoop
,
S. P.
Balaji
,
T. J. H.
Vlugt
, and
D.
Dubbeldam
,
J. Chem. Theory Comput.
10
,
942
(
2014
).
22.
F. A.
Escobedo
and
F. J.
Martinez-Veracoechea
,
J. Chem. Phys.
129
,
154107
(
2008
).
23.
F. A.
Escobedo
and
J. J.
dePablo
,
J. Chem. Phys.
105
,
4391
(
1996
).
24.
A. Z.
Panagiotopoulos
,
Int. J. Thermophys.
10
,
447
(
1989
).
25.
M.
Mehta
and
D. A.
Kofke
,
Mol. Phys.
79
,
39
(
1993
).
26.
D. A.
Kofke
and
E. D.
Glandt
,
Mol. Phys.
64
,
1105
(
1988
).
27.
A. Z.
Panagiotopoulos
,
N.
Quirke
,
M.
Stapleton
, and
D. J.
Tildesley
,
Mol. Phys.
63
,
527
(
1988
).
28.
J. J.
Depablo
and
J. M.
Prausnitz
,
Fluid Phase Equilib.
53
,
177
(
1989
).
29.
J. T.
Kindt
,
J. Chem. Phys.
143
,
124109
(
2015
).
30.
R. L. C.
Vink
and
J.
Horbach
,
J. Chem. Phys.
121
,
3253
(
2004
).
31.
C. M.
Wijmans
,
B.
Smit
, and
R. D.
Groot
,
J. Chem. Phys.
114
,
7644
(
2001
).
32.
Z.
Guo
and
J. T.
Kindt
,
Mol. Simul.
44
,
300
(
2018
).
33.
M. G.
Martin
and
J. I.
Siepmann
,
J. Am. Chem. Soc.
119
,
8921
(
1997
).
34.
B.
Chen
,
J. I.
Siepmann
, and
M. L.
Klein
,
J. Am. Chem. Soc.
124
,
12232
(
2002
).
35.
C. D.
Wick
,
J. I.
Siepmann
, and
D. N.
Theodorou
,
J. Am. Chem. Soc.
127
,
12338
(
2005
).
36.
P.
Bai
and
J. I.
Siepmann
,
J. Chem. Theory Comput.
13
,
431
(
2017
).
37.
J. M.
Stubbs
,
J. J.
Potoff
, and
J. I.
Siepmann
,
J. Phys. Chem. B
108
,
17596
(
2004
).
38.
39.
L.
Martinez
,
R.
Andrade
,
E. G.
Birgin
, and
J. M.
Martinez
,
J. Comput. Chem.
30
,
2157
(
2009
).
40.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
,
J. Mol. Graphics Modell.
14
,
33
(
1996
).
41.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
(
Academic Press
,
2002
).
42.
A.
Torres-Knoop
,
N. C.
Burtch
,
A.
Poursaeidesfahani
,
S. P.
Balaji
,
R.
Kools
,
F. X.
Smit
,
K. S.
Walton
,
T. J. H.
Vlugt
, and
D.
Dubbeldam
,
J. Phys. Chem. C
120
,
9148
(
2016
).
43.
W.
Shi
and
E. J.
Maginn
,
J. Comput. Chem.
29
,
2520
(
2008
).
44.
T. C.
Chu
,
R. J. J.
Chen
,
P. S.
Chappelear
, and
R.
Kobayashi
,
J. Chem. Eng. Data
21
,
41
(
1976
).
45.
I.
Wichterle
and
R.
Kobayashi
,
J. Chem. Eng. Data
17
,
9
(
1972
).
46.
I.
Wichterle
and
R.
Kobayashi
,
J. Chem. Eng. Data
17
,
4
(
1972
).
47.
J. M.
Beaudoin
and
J. P.
Kohn
,
J. Chem. Eng. Data
12
,
189
(
1967
).
48.
R. S.
Poston
and
J. J.
Mcketta
,
J. Chem. Eng. Data
11
,
362
(
1966
).
49.
H. L.
Chang
,
L. J.
Hurt
, and
R.
Kobayashi
,
AIChE J.
12
,
1212
(
1966
).
50.
J. P.
Kohn
and
W. F.
Bradish
,
J. Chem. Eng. Data
9
,
5
(
1964
).
51.
D. G.
Elliot
,
R. J. J.
Chen
,
P. S.
Chappelear
, and
R.
Kobayashi
,
J. Chem. Eng. Data
19
,
71
(
1974
).
52.
J. J.
Potoff
and
D. A.
Bernard-Brunel
,
J. Phys. Chem. B
113
,
14725
(
2009
).
53.
J. J.
Potoff
,
J. R.
Errington
, and
A. Z.
Panagiotopoulos
,
Mol. Phys.
97
,
1073
(
1999
).
54.
B.
Chen
and
J. I.
Siepmann
,
J. Phys. Chem. B
103
,
5370
(
1999
).
55.
C.
McCabe
,
A.
Galindo
,
A.
Gil-Villegas
, and
G.
Jackson
,
J. Phys. Chem. B
102
,
8060
(
1998
).
56.
S.
Aparicio
,
J. Supercrit. Fluids
46
,
10
(
2008
).
57.
J. D.
Haley
and
C.
McCabe
,
Fluid Phase Equilib.
440
,
111
(
2017
).
58.
J. H.
Simons
and
J. W.
Mausteller
,
J. Chem. Phys.
20
,
1516
(
1952
).
59.
L.
Zhang
and
J. I.
Siepmann
,
J. Phys. Chem. B
109
,
2911
(
2005
).
60.
H. A.
Lorentz
,
Ann. Phys.
248
,
127
(
1881
).
61.
D.
Berthelot
,
C. R. Hebd. Seanc. Acad. Sci. (Paris)
126
,
1703
(
1898
).
62.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
63.
G. C.
Boulougouris
,
I. G.
Economou
, and
D. N.
Theodorou
,
J. Phys. Chem. B
102
,
1029
(
1998
).
64.
E. W.
Lemmon
,
M. O.
McLinden
, and
D. G.
Friend
,
NIST Chemistry WebBook, NIST Standard Reference Database Number 69
, edited by
P. J.
Linstrom
and
W. G.
Mallard
(
National Institute of Standards and Technology
,
Gaithersburg MD
,
2000
), http://webbook.nist.gov.
65.
M. G.
Martin
and
A. L.
Frischknecht
,
Mol. Phys.
104
,
2439
(
2006
).
66.
A.
Sepehri
,
T. D.
Loeffler
, and
B.
Chen
,
J. Chem. Theory Comput.
13
,
1577
(
2017
).
67.
B. D.
Smith
and
R.
Srivastava
,
Thermodynamic Data for Pure Compounds: Part A Hydrocarbons and Ketones
(
Elsevier
,
Amsterdam
,
1986
).

Supplementary Material