Continuous fractional component Monte Carlo simulations of high-density adsorption in metal–organic frameworks

The continuous fractional component Monte Carlo method, which was designed to overcome difficulties with insertions and deletions of molecules, is modified to include configurational bias Monte Carlo methods and is further extended to binary systems. The modified method is shown to correctly predict adsorption of Ar in silicalite, Xe and Kr in HKUST-1, and enantiomers in a homochiral metal–organic framework. The modified method is also found to be approximately an order of magnitude more efficient in inserting and deleting molecules than traditional configurational bias grand canonical Monte Carlo simulations in dense systems.


Introduction
Metal -organic frameworks (MOFs) are a novel class of porous crystalline materials, built from metal corners or 'nodes' and organic 'linkers'. The great variety of nodes and linkers, and the large number of resulting possible combinations allow for tailoring and design of these materials towards particular applications. [1 -6] Molecular modelling has played an important role in the development, design, screening and understanding of MOFs for various applications such as methane storage, hydrogen storage, Xe/Kr separations, CO 2 capture and paraffin -olefin separations. [7 -13] Most of the focus has been on gas phase applications, but MOFs also have potential for liquid phase applications such as enantioselective separations, [14 -21] separation of industrial liquid mixtures [22 -27] and removal of organosulphur compounds from fuel. [28] Molecular modelling studies of liquid phase adsorption in MOFs are scarce. [17,18,20,21,29] Grand canonical Monte Carlo (GCMC) is the simulation method of choice to assess MOF performance for gas storage and gas separation, and GCMC is, in principle, applicable to liquid phase adsorption. The accuracy, precision and overall reliability of GCMC simulations depend heavily on the ability to insert and delete adsorbate molecules. When the systems of interest are too dense, as in liquid phase adsorption, insertion and deletion moves are rarely accepted, which leads to poor sampling of the loading and composition of the adsorbed phase.
Several biasing techniques have been developed to improve the acceptance rates for insertions and deletions in dense systems. Cavity biasing, [30] energy biasing, [31] configurational bias, [32 -34] parallel tempering [35 -37] and multi-canonical methods [38,39] are some of the techniques introduced to aid in the simulations of dense systems. Desgranges and Delhommelle have recently shown that an expanded Wang -Landau simulation method can be used to accurately calculate the entire adsorption isotherm of dense systems or large molecules by including the use of a fractional molecule. [40,41] Of particular interest to this study is further development of Monte Carlo algorithms using such fractional molecules. The continuous fractional component Monte Carlo (CFC MC) move, developed by Shi and Maginn, [42,43] is a method for the efficient insertion and deletion of molecules that can be applied in the grand canonical and osmotic ensembles. The CFC move also draws from various other methods. [44 -52] In this study, we extend the CFC MC move for use in mixtures, focusing on adsorption in MOFs. We first validate our extension of the move to binary mixtures of Xe and Kr in the MOF HKUST-1. Finally, we study the adsorption of racemic mixtures of R-and S-dimethylallene (DMA) and racemic mixtures of R-and S-dimethylcyclobutane (DMB) in the homochiral MOF, Cd-BINOL. [53] Using CFC MC, we show that it is possible to simulate adsorption of (R,S)-DMA and (R,S)-DMB in Cd-BINOL even at high adsorbed phase densities.

Methodology
The grand partition function for an adsorbed system with fixed temperature T, chemical potential m and volume V is where b ¼ 1/k B T, N is the number of molecules in the system, Z(N,V,T) is the configurational integral and q contains kinetic energy contributions from all the atoms in the system and the Jacobian of the transformation from atomic Cartesian coordinates to internal coordinates. CFC MC expands on the original grand partition function (Equation (1)) by including a tagged particle, or 'fractional' molecule, with scaled interactions using a parameter l. l is a pseudo-continuous coupling parameter that is used to scale interactions between the fractional molecule and real molecules by modifying the standard 12-6 Lennard Jones potential into the following functional form: ; ð2Þ where f f ðr if ; lÞ is the potential energy between real atom i and fractional atom f, r if is the distance between atom i and fractional atom f, and s and 1 are the Lennard-Jones parameters.
For a system with more than one component, the system can be expanded with as many tagged molecules as components following the approach of Escobedo and de Pablo to give the extended grand partition function [44]: where the total number of molecules . . . ; j t ; V; TÞ is again the configurational integral, q now contains kinetic energy contributions from all atoms (real and fractional) in the system and the Jacobian of the transformation from atomic Cartesian coordinates to internal coordinates and t is the number of fractional molecules which equals the number of components. Each fractional molecule can exist in a large number of states m i with each state j i having an expanded system variable h j i associated with it. h j i is the bias factor associated with the particular value of l i at state j i . In systems with only one component and one fractional molecule, Equation (3) In order to verify new stochastic moves, we need to show that the new move obeys microscopic reversibility by imposing the strict condition: where P mn is the one-step transition or acceptance probability of going from state m to state n and r m and r n are the ensemble probability distribution functions for states m and n, respectively. By working with fugacities (f ) instead of chemical potentials we can show that the probability distribution function for some state m having N m real molecules and t number of components and fractional molecules is [42,50] where the total number of molecules in state m is f m is the total potential energy of state m, V t is the volume associated with the generalised coordinates of a single molecule of species t and Z IG t is the ideal gas configurational integral of a single molecule of species t.
Microscopic reversibility in Equation (5) is achieved by using the following acceptance rule: where a mn is the probability of attempting a move from state m to state n. For systems that contain more than one fractional molecule, we use a separate value of l for each fractional molecule. Going through a similar derivation as that used by Shi and Maginn [42] for the single fractional molecule probabilities, we arrive at a similar result for the changes in l for mixtures with t fractional molecules. In the original CFC MC method, fractional molecules are always present in the system, and attempts to change l are made as Monte Carlo moves. During a change in l, if l becomes greater than 1 then the fractional molecule is converted to a real molecule and a new fractional molecule is inserted randomly into the system with l new being equal to l 2 1. The new state n contains N n ¼ N m þ 1 real molecules and t fractional molecules. The move is accepted with a probability similar to that of the grand canonical ensemble insertion move by substituting Equation (6) into Equation (7): where Df mn is the change in total potential energy when going from state m to state n and f , V and Z IG correspond to the species being inserted. If l becomes less than 0, then the fractional molecule is removed and one of the remaining real molecules of the same component is randomly chosen to become the new fractional molecule with l new equal to l þ 1. The move from state m to state n is accepted with a probability similar to that of the grand canonical ensemble deletion move: Again f , V and Z IG correspond to the species being deleted. For moves in which changes in l do not either exceed 1 or fall below 0, the acceptance probability is where, for this move, the a mn and a nm attempt probabilities are equal and cancel out and are not explicitly shown. Torres-Knoop et al. [52] were able to demonstrate that in dense systems, the acceptance rate of certain moves can be increased by allowing fractional molecules to undergo configurational bias Monte Carlo (CBMC) moves. For example, a fractional molecule could be 'reinserted' by removing it and regrowing it using configurational bias. In addition, for changes in l where l goes above 1, the new fractional molecule could be grown into a low energy position using configurational bias in a combined CFC CB MC move. The approach of Torres-Knoop et al. also eliminates the need for generation of ideal gas configurations and allows the fractional molecule's configuration to adopt configurations that could deviate from its ideal gas configuration. [52] In this study, we attempted CFC CB insertion and reinsertion moves for the fractional molecules as described by Torres-Knoop et al. [52]. A reinsertion move removes a molecule and attempts to reinsert it at a random location in the system. Reinsertion moves differ from standard translation moves in that standard translation moves are performed to achieve a 50% acceptance rate by adjusting the maximum distance a molecule can translate from its original position. For the CFC insertion moves (i.e. changes in l where l goes above 1), we can use derivations from Torres-Knoop et al. [52] and from other work, [54,55] to include CBMC probabilities in the ratio of insertion attempt probabilities, in the ratio of deletion attempt probabilities, and for the ideal gas Rosenbluth weight where f ext is the intermolecular energy of the growing molecule, C is a normalisation constant, [55] v is the Rosenbluth weight and v IG is the external average ideal gas Rosenbluth weight, where the only interactions that contribute to v IG are the intramolecular non-bonded interactions. CBMC is often applied to long, flexible chain molecules, but as described in the following paragraph, it is also useful for improving insertion acceptance rates for rigid molecules. In this work, we have restricted ourselves to rigid molecules. Substituting Equation (11) into Equation (8), Equation (12) into Equation (9) and noting that Equation (13) reduces to Z IG =ðVCÞ ¼ 1 for monatomic species and rigid molecules: As above, f corresponds to the species being inserted or deleted.
In this way CBMC techniques can be implemented in tandem with CFC MC moves. For rigid molecules, only three atoms need to be placed to define the location and orientation of the molecule. For molecules that are newly inserted into the system, the so-called high loading firstatom sampling [55] is performed. For the first atom, n trial insertions are attempted at random locations, with the total potential energy of the atom calculated at each location. A location r i , that is weighted towards lower energy positions, is then chosen with the following probability: and Rosenbluth value: The placement of the second and third atoms uses the methods of sphere sampling and disc sampling, respectively. [55] The probabilities and corresponding Rosenbluth weights (v sphere and v disk ) for these moves can be found in the work of Macedonia and Maginn. [55] The new molecule's Rosenbluth weight is the product of the previously defined Rosenbluth weights: which can be substituted into Equations (14) and (15). Sometimes l can become 'stuck', mimicking similar problems that the use of the CFC GCMC method is trying to solve. This can occur when the system's energetics prefer that the fractional molecule exists in a small range of l values, and subsequent Monte Carlo moves that attempt to change l to values outside this small range are difficult to accept. In order to overcome this problem, we can use a Wang -Landau weighting factor [56,57] for h to make the probability of states nearly equal for all values of l. This is accomplished by dividing the l range into 10 bins and setting the initial biasing factors to 0. During equilibration of the system, the values of l are recorded and binned accordingly. The bin's weighting factor h is then adjusted to h ¼ h 2 lnðvÞ, where v is the scaling parameter and initially lnðvÞ is unity. Every 10,000 CFC MC steps during the equilibration cycles, the distribution of l was checked and if the probability of each l bin was greater than 1%, v was modified to v ¼ ffiffiffi v p .

Simulation details
Before the simulation, in order to speed up calculations, framework -adsorbate potentials can be tabulated on three-dimensional grids throughout the unit cell. During the simulation, the potential energies between the framework and adsorbate atoms are determined from the pretabulated information using a three-dimensional cubic hermite polynomial interpolation algorithm. [58] For this work, to allow fractional molecules to be used with pretabulated grids, l was considered as unity for fractionalframework intermolecular interactions. In addition, to prevent fractional molecule overlap with other fractional molecules, l was considered as unity for all fractionalfractional intermolecular interactions. All Lennard-Jones interactions were scaled using Equation (2). None of the systems considered in this work included Coulombic interactions. All simulations using CFC MC moves used the Wang -Landau weighting (described in the previous section). Note that molecules could attempt more than one type of MC move and that all molecules were picked from their corresponding MC moves with equal probability. All simulations were carried out using our in-house code RASPA. [59]

Argon in silicalite
The CFC CBMC move was first implemented for simulations of argon in the zeolite silicalite. Three types of simulations were performed: a CB GCMC simulation with no fractional atoms, a CB GCMC simulation using only CFC moves to add and remove atoms from the system and a CB GCMC simulation using CFC moves but also using standard insertion and deletion moves for real atoms.
In addition to the moves specified, all simulations used translation and reinsertion moves. A pre-tabulated grid for Ar -O interactions (0.1 Å spacing), a cut-off of 13.0 Å and 2 £ 2 £ 2 unit cells were used. All silicalite atoms were held fixed in their crystallographic positions. The simulations used 1,000,000 equilibration and 1,000,000 production cycles. A cycle consists of the minimum of 20 and N Monte Carlo steps, where N is the number of argon atoms in the system at the beginning of the cycle. This work used the Lennard-Jones parameters of Snurr et al. [31] where only the Ar -O and Ar-Ar interactions were considered. The parameters can be found in Table S1 in the supplemental data (SD).

Xenon and krypton in HKUST-1
To test the CFC MC move for adsorption simulations of binary mixtures, pure Xe, Kr and an 80/20 Kr/Xe mixture were simulated in the MOF HKUST-1 [60] at low and high densities at 273 K. Insertion and deletion moves were only used in the traditional CB GCMC simulations. In binary mixture simulations, an identity change move was allowed, where a species of component A could attempt to convert to component B or vice versa. All systems were allowed translation and reinsertion moves. For the binary systems containing fractional atoms, four different types of CFC simulations were tested. In the first case, CFC moves were attempted only for Kr. In the second case, CFC moves were attempted only for Xe. In the third case, one fractional atom was present and it was allowed to undergo identity change moves. Depending on its identity, CFC moves were then attempted for either Xe or Kr. In the last case, two fractional atoms were present, one Xe and one Kr. Both the traditional CB GCMC simulation and the CB GCMC simulations using CFC moves were performed using 100,000 equilibration and 100,000 production cycles. We used the Lennard-Jones parameters of Ryan and co-workers [61], which can be found in Table S2 in the SD. A pre-tabulated grid for framework -adsorbate interactions (0.1 Å spacing), a cut-off of 12.0 Å and a 1 £ 1 £ 1 unit cell were used. As in silicalite, all framework atoms were fixed in their crystallographic positions.

(R, S)-DMA and (R, S)-DMB in Cd-BINOL
Finally, we performed CFC CB GCMC simulations for molecules. Specifically, racemic mixtures of (R,S)-DMA and racemic mixtures of (R,S)-DMB were simulated in the MOF Cd-BINOL [53], following the work of Bao et al. [17]. All simulations used translation, rotation and reinsertion moves. Insertion and deletion moves were exclusive to the traditional CB GCMC simulations. For the CFC CB GCMC binary simulations, four simulations were performed as described above for Xe/Kr. In this case, the two components are R-DMA and S-DMA or R-DMB and S-DMB. Both the traditional CB GCMC simulation and the CFC CB GCMC simulations were performed using 100,000 equilibration and 100,000 production cycles. The Lennard-Jones parameters were taken from the work by Bao and co-workers [17] and are given in Table S3 in the SI. A pre-tabulated grid for framework -adsorbate interactions (0.1 Å spacing), a cut-off of 12.8 Å and 2 £ 2 £ 1 unit cells were used. All framework atoms were fixed in their crystallographic positions.

Argon in Silicalite
A series of CFC CB GCMC and CB GCMC simulations were conducted for Ar in silicalite at 77 K and various pressures (0.2665, 1.01 and 101.3 Pa). The results were compared against those previously reported by Snurr and coworkers [31] and are found in Table 1. The results from Snurr et al. agree, as expected, with our traditional CB GCMC simulations and the CFC CB GCMC simulations at all loadings. We also found good agreement for the average potential energy per atom kf=Nl in all cases. The good agreement for both the loading and energies indicates that the CFC algorithm has been correctly implemented in our software. Although traditional CB GCMC simulations would suffice for simulating Ar in silicalite at these state points, these simulations were used to test the CFC CB GCMC method under conditions where traditional CB GCMC will give reliable results to compare against. These results also demonstrate that CFC CB GCMC can accurately predict adsorption of gas species in an adsorbent.

Xenon and krypton in HKUST-1
To extend the work to mixtures, we performed both CFC CB GCMC and CB GCMC simulations for Xe and Kr in HKUST-1 at 273 K and at three different pressures (0.1, 1 and 3 MPa). We compared our results from traditional CB GCMC and CFC CB GCMC with two fractional atoms against those of Ryan et al. [61] and found good agreement for all cases, as shown in Figure 1. Single component isotherm results can be found in Figure S1 in the SI and show similar agreement.
To further test our CFC CB GCMC method, we examined different ways of treating the fractional atom(s) in binary systems, using four different options: one Kr fractional atom only, one Xe fractional atom only, one fractional atom that was allowed to change identity during the simulation and two fractional atoms (one Kr and one Xe). The tabulated results for all four options can be found in Table 2, where we found good agreement for all cases. This indicates that the CFC CB GCMC method can be correctly implemented in different ways for binary mixtures.  Figure 2 shows the total energy of the system per real atom as a function of pressure. We notice that all systems simulated with one or more fractional atoms have systematically lower energies, which we attribute to the fractional atom(s) in the system. This systematic error is most noticeable at low loadings, where frameworkfractional and real -fractional interactions (which are included in the total system energy) account for a larger amount of the total energy. At higher loadings, where the contribution from the fractional atom is a smaller portion of the total energy, this discrepancy decreases. This supports the idea that we can relax our physical picture of the system and allow configurations with fractional atoms to be valid representations of the system at higher loadings, as suggested by Escobedo and de Pablo [44].

(R, S)-DMA and (R, S)-DMB in Cd-BINOL
After verifying that CFC CB GCMC worked for both single component and binary monatomic species, we applied the method to simulate racemic mixtures of (R,S)-DMA in Cd-BINOL. The CFC CB GCMC method was first compared against the traditional CB GCMC method at seven different pressures (50, 150, 300, 500, 1000, 2000 and 4000 Pa), as shown in Figure 3. We find that all three methods (traditional CB GCMC, CFC CB GCMC using one fractional molecule and CFC CB GCMC using two fractional molecules) are in good agreement over the entire range. Similar results were obtained for racemic mixtures of (R,S)-DMB. The results indicate that the CFC CB GCMC method is working correctly and can now be used for systems with higher densities where traditional CB GCMC simulations are inefficient or fail.
Results for (R,S)-DMA in Cd-BINOL at 2 MPa and varying bulk phase composition are shown in Figure 4, where we have plotted the loading of R and S molecules and the change in enantiomeric excess (Dee ¼ ee MOF 2 ee bulk ) against the bulk phase mole fraction of R-DMA. The first thing we notice is that the number of adsorbed    Table S5 (SD). molecules of R-and S-DMA in the 0.50 R-DMA mole fraction mixture at high densities is not equal, and is actually equal at a lower mole fraction of R-DMA, around 45%. This can be seen where R-DMA crosses S-DMA in Figure 4. This lower-than-expected mole fraction indicates that even at high densities Cd-BINOL has a slight preference to adsorb R-DMA over S-DMA, as shown in Figure 4. We further find that the enantiomeric excess (ee) decreases when going from gas phase densities (ee ,60% -20%) [17] to liquid phase densities (ee ,10%). The fact that there is a high ee at low densities and a low ee at high densities further supports the finding of Bao and coworkers [17] that Cd-BINOL has specific enantioselective sites and that the majority of these sites become occupied at low pressures. Upon increasing the mole fraction of R-DMA, we find that Dee decreases, suggesting that as the system becomes saturated with R-DMA, the sites in Cd-BINOL have less effect on the ee. Figure 4 also shows that the loadings for the CFC method and traditional CB GCMC agree over the entire range of mole fractions of R-DMA. However, we found that the CFC CB GCMC simulations are more efficient than traditional CB GCMC simulations for both the system containing (R,S)-DMA and the system containing (R,S)-DMB, as shown in Figure 5, where efficiency is defined as the number of accepted insertion and deletion moves of real molecules per total simulation time, or, how many insertions or deletions are occurring per second of simulation time. We find that the CFC CB GCMC method is approximately 6 to 10 times more efficient than the traditional CB GCMC method for (R,S)-DMA. For the (R, S)-DMB system, the CFC CB GCMC method was approximately an order of magnitude more efficient than the CB GCMC method. We expect this ratio may be even higher for other systems, especially those where traditional CB GCMC simply fails.

Conclusion
The CFC CB GCMC method was extended to mixtures and applied to simulate adsorption in several adsorbents. We first validated the method for single component Ar in silicalite and binary mixtures of Xe and Kr in the MOF HKUST-1. For the binary system, we tested three different approaches: (1) a system where only one fractional atom was present and the fractional atom maintained the same identity throughout the simulation, (2) a system with only one fractional atom present and the fractional atom was allowed to change its identity and (3) a system with two fractional atoms, one fractional atom for each of the binary components present in the system. In all cases, we found good agreement in loadings between the CFC CB GCMC method and the traditional CB GCMC method for both the single component and binary mixture systems. However, we found that the CFC CB GCMC method systematically predicted slightly lower energies than traditional CB GCMC methods. This small discrepancy in energies decreased as the number of molecules in the system increased. For our high-density simulations, we expect the error in energy from the fractional molecule to be around 2%. We further modified the CFC CB GCMC method to incorporate configuration biasing for use with molecules. The CFC CB GCMC method was tested using a racemic mixture of (R,S)-DMA and a racemic mixture of (R,S)-DMB in the MOF Figure 4. Adsorption of (R,S)-DMA in Cd-BINOL at high density (P ¼ 2 MPa, T ¼ 300 K). Results from traditional CB GCMC and CFC CB GCMC, using two fractional molecules, are compared at varying mole fractions of R-DMA in the bulk phase. Change in enantiomeric excess (Dee) from the CFC CB GCMC simulation is represented by the dashed line. Tabulated values can be found in Table S6 (SD). Cd-BINOL, and we found similar results for a variety of gas phase pressures compared with the traditional CB GCMC method.
After the CFC CB GCMC method was validated for mixtures, it was used to simulate a high-density system of (R,S)-DMA and a high-density system of (R,S)-DMB in Cd-BINOL, where we found a sharp decrease in ee when moving from gas phase adsorption to liquid phase adsorption. We also found that CFC CB GCMC was ,8 times more efficient at inserting and deleting molecules than the traditional CB GCMC.

Disclosure statement
No potential conflict of interest was reported by the authors.

Supplemental data
Supplemental data for this article can be accessed here