Using Monte Carlo Simulation to Compute LiquidâVapor Saturation

Article pubs.acs.org/JPCB

Using Monte Carlo Simulation to Compute Liquid−Vapor Saturation Properties of Ionic Liquids Kaustubh S. Rane and Jeffrey R. Errington* Department of Chemical and Biological Engineering, University at Buffalo, The State University of New York, Buffalo, New York 14260-4200, United States S Supporting Information *

ABSTRACT: We discuss Monte Carlo (MC) simulation methods for calculating liquid− vapor saturation properties of ionic liquids. We first describe how various simulation tools, including reservoir grand canonical MC, growth-expanded ensemble MC, distance-biasing, and aggregation-volume-biasing, are used to address challenges commonly encountered in simulating realistic models of ionic liquids. We then indicate how these techniques are combined with histogram-based schemes for determining saturation properties. Both direct methods, which enable one to locate saturation points at a given temperature, and temperature expanded ensemble methods, which provide a means to trace saturation lines to low temperature, are discussed. We study the liquid−vapor phase behavior of the restricted primitive model (RPM) and a realistic model for 1,3-dimethylimidazolium tetrafluoroborate ([C1mim][BF4]). Results are presented to show the dependence of saturation properties of the RPM and [C1mim][BF4] on the size of the simulation box and the boundary condition used for the Ewald summation. For [C1mim][BF4] we also demonstrate the ability of our strategy to sample ion clusters that form in the vapor phase. Finally, we provide the liquid− vapor saturation properties of these models over a wide range of temperature. Overall, we observe that the choice of system size and boundary condition have a non-negligible effect on the calculated properties, especially at high temperature. Also, we find that the combination of grand canonical MC simulation and isothermal−isobaric temperature expanded ensemble MC simulation provides a computationally efficient means to calculate liquid−vapor saturation properties of ionic liquids.

1. INTRODUCTION Room temperature ionic liquids (RTILs) are considered as potential next generation industrial solvents due to their unique properties. Knowledge of liquid−vapor coexistence properties of these fluids is of general interest. Molecular simulation can play an important role in the rational design of new ionic liquids and in studying the molecular aspects of poorly understood phenomena. However, their application to the study of liquid−vapor phase behavior has proven to be challenging. In this work, we show how Monte Carlo simulation strategies recently developed within our group1,2 are used to calculate the liquid−vapor coexistence properties of ionic liquids. We present results for the restricted primitive model (RPM) and a realistic model for 1,3-dimethylimidazolium tetrafluoroborate ([C1mim][BF4]) over a broad range of temperature. Molecular simulations have greatly contributed to the study of RTILs. For example, they facilitate the calculation of properties at temperatures higher than their thermal degradation temperature. This enables one to obtain estimates for critical parameters, which are important for the application of corresponding states theories.3−5 Recently, Rai and Maginn used Gibbs ensemble Monte Carlo simulations employing configurational-bias sampling to calculate liquid−vapor saturation properties for certain ionic liquids.3,4 These authors considered the homologous series of imidazolium based RTILs © XXXX American Chemical Society

with different alkyl chain lengths associated with the cation. The properties that were obtained included critical parameters, vapor pressures, and enthalpies of vaporization. They highlighted the effect of temperature and alkyl chain length on these quantities, and also compared the observed trends with experimental data. These studies demonstrated that molecular simulations performed using available force-fields for RTILs provide reasonable descriptions of how the nature of interionic interactions influences phase coexistence properties. Experiments have also been used to measure the liquid− vapor saturation properties of RTILs.6−8 Rocha et al. studied the vapor pressure and enthalpy and entropy of vaporization for the 1-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide series.6 However, the range of temperatures investigated by experimental and Monte Carlo simulation studies rarely overlap. Monte Carlo simulations generally encounter sampling difficulties at low to moderate temperatures, which have limited their applicability to relatively high temperatures. Experiments, on the other hand, are performed at moderate temperatures (400−600 °C) due to the chemical decomposition of RTILs at higher temperatures and technological difficulties in measuring vapor phase properties at room temperature.6,8 We show here Received: April 29, 2013 Revised: June 3, 2013

A

dx.doi.org/10.1021/jp404207x | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

system sizes for future simulation studies of ionic liquids. It has been observed previously that details of the boundary condition used in the Ewald summation influences the system size dependence of critical parameters of RPM calculated from GC simulations.20 In this work, we therefore employ both the tinfoil and vacuum boundary conditions to study finite-size effects on saturation properties. The computation of phase coexistence properties of simple models of ionic liquids enjoys a rich history.5,12,20−24 The simplest model of ionic liquids is the restricted primitive model (RPM), which consists of oppositely charged hard spheres with equal diameters. The liquid−vapor coexistence densities, pressures, and critical parameters of RPM have been reported by several groups12,20,24−27 and sufficient data are available in the literature to perform benchmark validation tests for new simulation strategies applied to ionic systems. We present results for the phase coexistence properties of RPM over a wide range of temperatures obtained from grand canonical temperature expanded ensemble simulation2 and compare these with those from earlier studies. In comparison to RPM, very few molecular simulation studies have examined the liquid−vapor coexistence properties of realistic models for ionic liquids.3,4 In this paper, we show the results for a united-atom model of [C1mim][BF4], which was studied by Rai and Maginn.3 We calculate the saturated densities and vapor pressures for this fluid at temperatures ranging from 300 to 1200 K using a combination of direct grand canonical and isothermal−isobaric temperature expanded ensemble simulation. The paper is organized as follows. In the following section we describe the simulation methods used to sample systems of interest and the techniques that are used to compute phase coexistence properties. Here, we place emphasis on describing simulation tools important for ionic systems. Within Section 3 we discuss the molecular models employed and provide details for the simulations performed within this study. In Section 4 the results obtained are presented and compared with those in the literature. The paper is concluded in Section 5.

that Monte Carlo simulation based strategies introduced recently1,2 provide a means to compute saturation properties of model RTILs at experimentally accessible temperatures. Our approach uses a combination of direct and temperature expanded ensemble simulations.2 Direct grand canonical (GC) or isothermal−isobaric ensemble (NPT) Monte Carlo simulations are performed to calculate saturation properties at select temperatures. These are followed by either grand canonical or isothermal−isobaric temperature expanded ensemble (TEE) simulations that separately track the Gibbs free energy or grand potential, respectively, of the coexisting phases along the saturation curve. Histogram reweighting schemes allow us to locate saturation conditions where Gibbs free energies (grand potentials) of the two phases are equal. This strategy helps in overcoming difficulties often encountered while sampling densities intermediate between those characteristic of the liquid and vapor phases. To apply this strategy to complex molecules, it is necessary to ensure proper sampling of the intramolecular configuration space. To address this challenge we use simulation tools like reservoir grand canonical Monte Carlo,9 single molecule hybrid Monte Carlo (HMCOne), and conventional hybrid Monte Carlo (HMC) moves.10 The performance of GC simulations greatly depends on the acceptance of insertion and deletion moves. We employ techniques like growth expanded ensemble11 moves to address difficulties encountered in inserting and deleting complex molecules from the dense liquid phase. The approach noted above has been successfully applied to study liquid−vapor phase equilibria of various complex molecules, such as cyclohexane, pyrene, and squalane.2 Simulations of ionic liquids encounter additional challenges like strong association of ions at distances close to contact and formation of aggregates at low densities and temperatures. To address these issues, Orkoulas and Panagiotopoulos introduced distance-biasing particle insertion/deletion and cluster moves.12 Recently, these moves have been used in Gibb’s ensemble simulations of imidazolium based ionic liquids.3,4 In this work, we use distance−biasing moves with a slightly different biasing strategy to account for the flexibility of the molecules. Additionally, these moves are combined with reservoir grand canonical MC and growth expanded ensemble strategies to facilitate the insertion and deletion of ions in the liquid phase.2 To validate our simulation tools we compare results for [C1mim][BF4], propanol, and octane with and without distance-biasing moves. Since Earle et al.13 demonstrated that ionic liquids can be distilled, great interest has been expressed in understanding the nature of their vapor phase.14,15 For example, it has been observed that an increase in temperature favors the formation of larger aggregates and ion-pairs are predominant at low temperatures. To study the vapor phase over a wide range of temperature, it is necessary that the simulation tools allow proper sampling of ion clusters of variable size. In our simulations we employ an extension of aggregation-volume-bias Monte Carlo moves (AVBMC3) proposed by Chen and Siepmann16 to allow formation and breakage of ion clusters. The limitation of computational resources places an upper limit on the size of the system used to perform simulations. It is well-known that liquid−vapor coexistence properties can show non-negligible system size dependence, especially at temperatures near the critical point.17−20 Here, we provide results from finite-size studies performed with the RPM and [C1mim][BF4]. Such studies help us to identify appropriate

2. SIMULATION METHODS We first outline the molecular simulation tools that are used within this study. The features of many of these techniques are detailed within an earlier report.2 We highlight those features that are most relevant to the simulation of ionic liquids. The main challenge in performing Monte Carlo simulations of these fluids is the sampling of low density and low temperature states that are characterized by the presence of ion pairs or clusters. In addition to strongly associating character, RTILs have complex structures with intraionic degrees of freedom which provide additional challenges when performing molecular transfers in grand canonical simulations. In the following paragraphs we discuss the strategies used to address these issues. Reservoir Grand Canonical MC. Within a grand canonical (GC) simulation, individual ions are exchanged between a system of interest at a given temperature and activity and a reservoir of ions at the same temperature. In this work, we use separate ideal gas reservoirs for cations and anions. Details regarding maintenance of the reservoirs are provided in an earlier report.2 Each reservoir contains a collection (say 500) of ions with unique intraionic configurations. The probability res pres i (ri ) of selecting a reservoir ion of type i with configuration res ri is given by B


The Journal of Physical Chemistry B res

pires (r ires)

Article

− j using the distance-biasing scheme relative to the random deletion case is

res

e−βEi (r i = zires

)

(1)

Eres i

Pdel = Nc pjcounter (rijcom)

zres i

where is the intraionic interaction energy and is the ionic configurational partition function. Distance-Biasing Insertion/Deletion. To insert an ionpair we use a distance-biasing algorithm proposed by Orkoulas and Panagiotopoulos.12 In this technique, oppositely charged ions are preferentially inserted near each other with a probability that is generally an exponentially decreasing function of the interionic distance.12,28 In this work, the center-of-mass separation distance rcom between ion i and ij counterion j is selected with a probability density W (rijcom)

⎧ 1 −β ⎪ e = ⎨C ⎪0 ⎩

U (rijcom)

In this work, we consider systems consisting of monovalent ions, and therefore Nc = N, where N is the number of ion-pairs. Growth Expanded Ensemble. We transfer ions to and from the system via an M-stage procedure in which the Lennard-Jones and Coulomb interactions between the selected (tagged) ion-pair and the rest of the system scale with the stage number m, while the two tagged ions interact with each other via the full potential. Within the growth expanded ensemble, a given configuration consists of N complete molecules and at most one partial molecule, where each molecule is composed of an ion-pair consisting of an ion of type i and an ion of type j. The probability πk of observing a given microstate k at temperature T and activity ξ is

rmin ≤ rijcom ≤ rmax rijcom ≤ rmin , rijcom ≥ rmax

(2)

where ⟨U(rcom ij )⟩ is the Boltzmann-averaged interionic interaction energy averaged over a suitable number of intraionic configurations and relative orientations for the given ion types. The normalization constant C is given by C=

∫r

rmax

e−β U (r) dr

min

πk =

(3)

4πr 2 V

⎧ ⎫ ⎪ π 1 P ⎪ pacc = min⎨1, n res res del ⎬ ⎪ ⎪ ⎩ πo pi pj Pins ⎭

(4)

W (rijcom) W0(rijcom)

⎧ π P ⎫ pacc = min⎨1, n pires pjres ins ⎬ Pdel ⎭ ⎩ πo

(5)

To delete or displace an ion-pair, the first ion i is selected randomly and the counterion j is selected with probability pjcounter (rijcom) =

N

(9)

(10)

where Pins now refers to the bias associated with inserting the selected ion-pair into microstate n and Pdel refers to the bias for removing the selected ion-pair from microstate o. We now consider the case in which ion-pairs are transferred between the system and the reservoir in a multistep process (M > 1). The acceptance probability for inserting a tagged ion-pair is given by

W (rijcom) ∑k =c 1 W (rikcom)

(8)

where Pins refers to the bias associated with inserting the trial ion-pair into microstate o and Pdel refers to the bias for removing the trial ion-pair from microstate n. The corresponding transition probability for removing a complete ion-pair is

where V is the volume of the simulation cell. The bias Pins associated with inserting the ion-pair i − j using the distancebiasing scheme relative to the random insertion case is Pins =

N +Δ 2 1 ⎛⎜ 1 ⎞⎟ 2(N +Δ)⎛⎜ ξ ⎞⎟ V e − βE k ⎜ z resz res ⎟ Ξ ⎝ N! ⎠ ⎝ i j ⎠

where Ξ is the partition function of the growth expanded ensemble, β = 1/kBT is the inverse temperature (kB is Boltzmann’s constant), Ek is the total configurational energy, and Δ is a parameter that evaluates to 0 when the system consists of full molecules only and 1 when the system contains res βμ a tagged molecule. The activity is defined as ξ = qizres i qjzj e , where μ is the chemical potential of the ion-pair and qi represents the component of the molecular partition function of ion i stemming from integration over the momenta. During a growth expanded ensemble MC simulation one attempts to insert (remove) an ion-pair, increase (decrease) the stage number m while the molecule remains partially grown, and complete (start) the growth (destruction) of an ion-pair. We first consider the case in which complete ion-pairs (M = 1) are transferred between the system and the reservoir. Let us say that one of these moves takes the system from microstate o to microstate n. The acceptance probability for inserting a complete ion-pair consisting of ions of types i and j is given by

The distance rmin is related to the repulsive core for the two ions and is established by locating the center-of-mass separation distance at which ⟨U(rcom ij )⟩ becomes sufficiently large and positive, thereby rendering W negligible. The distance rmax is selected for computational efficiency and is always taken to be less than half of the box size. We note that the functional form 28 for W(rcom ij ) is similar to that employed by Hynninen et al. com Also, we find that the effective potential ⟨U(rij )⟩ does not vary significantly with temperature, and therefore we use a single version of this potential for a given ion pair, which is evaluated at a temperature slightly below the expected critical temperature. During the insertion of an ion-pair consisting of ion i and counterion j, the position of the center-of-mass and orientation of ion i are selected randomly and the separation distance rcom ij com is selected with frequency W(rcom ij )/W0(rij ) relative to the random insertion of ion j. Here, W0(r) is the probability density of randomly placing two ions with a center-of-mass separation distance r within the simulation cell.12 Provided r is less than half of the simulation box length W0(r ) =

(7)

(6)

where Nc is the total number of counterions available within the simulation cell. The bias Pdel associated with selecting ion-pair i C


The Journal of Physical Chemistry B pacc

⎧ ⎫ ⎪ πn 1 1 ⎪ ⎨ ⎬ = min 1, res res ⎪ ⎪ ⎩ πo pi pj Pins ⎭

Article

and disruption of another in a single move. Here, two ions of the same type, k1 and k2, are selected such that their neighborhoods (bonded regions) are nonoverlapping. An ion i, which is of the same type as k1 and k2, is now selected using one of the following schemes: (1) i is randomly selected from the ‘in’ region of k2 and transferred to a random location within the ‘in’ region of k1, (2) i is randomly selected from the ‘out’ region of k1 and is transferred to the ‘in’ region of k1, (3) i is randomly selected from the ‘in’ region of k1 and is transferred to the ‘in’ region of k2, or (4) i is randomly selected from the ‘in’ region of k1 and is transferred to the ‘out’ region of k1. Note that the ‘out’ region of k1 also includes the ‘in’ region of k2, but excludes ion k2. Each of the four moves outlined is selected with equal probability (this corresponds to the Pbias = 0.5 case within ref 16). Next, the ion-pair is completed by selecting a counterion j according to the distance-biasing approach, i.e., of all the available counterions, j is selected with probability pcounter . As before, we will consider a transition occurring from j microstate o to microstate n. The acceptance probability for move (1) is given by

(11)

and the corresponding transition probability for removing a tagged ion-pair is ⎧ π ⎫ pacc = min⎨1, n pires pjres Pins⎬ ⎭ ⎩ πo

(12)

The ion-pair selection bias Pdel does not play a role in the two expressions above due to the fact that the tagged ions are selected for removal with a probability of one. Once a tagged ion-pair has been added to the system, the tag remains until the tagged ion-pair is removed or the insertion process is completed (m = M). The acceptance probability for a stage change move wherein the molecule remains tagged is simply ⎧ π ⎫ pacc = min⎨1, n ⎬ ⎩ πo ⎭

(13)

The transition probability for completing the growth of a tagged ion-pair (m = M in microstate n) is ⎧ π ⎫ pacc = min⎨1, n Pdel ⎬ ⎩ πo ⎭

pacc

Nkin1

(16)

Nkin2

where and correspond to the number of ions of type i in the ‘in’ regions of ions k1 and k2, respectively. The terms pcounter j,o and pcounter represent the probabilities of selecting counterion j j,n in the o and n microstates, respectively. The acceptance probability for move (3) is of the same form as in eq 16 with the k1 and k2 labels switched. The move (2) is accepted with probability

(14)

and that of initiating destruction of a complete ion-pair (m = M in microstate o) is ⎧ π 1 ⎫ ⎬ pacc = min⎨1, n ⎩ πo Pdel ⎭

⎧ ⎫ pjcounter ⎪ ⎪ πn Nink 2 ,n ⎬ = min⎨1, counter ⎪ k1 ⎪ ⎩ πo (Nin + 1) pj , o ⎭

(15)

counter ⎫ ⎧ k1 ⎪ Nout π Vin pj , n ⎪ ⎬ pacc = min⎨1, n counter ⎪ k1 ⎪ ⎩ πo (Nin + 1) Vout pj , o ⎭

Following the rationale provided above, the ion-pair insertion bias Pins does not play a role in eqs 14 and 15 because the tagged ions are selected for growth completion with a probability of one. Within Appendix A we provide detailed working expressions for eqs 9 to 15 that we expect to be useful in practice. Finally, we note that we schedule the number of growth stages M to increase with the total molecule number N for direct GC simulations or with the inverse temperature β for GC-TEE simulations. AVBMC3 and AVBMC2. The strongly associating nature of ionic liquids is responsible for difficulties in properly sampling relevant microstates using conventional MC moves. These sampling difficulties are particularly troublesome at low temperatures wherein the thermal energy of molecules is insufficient to disrupt aggregates. It is therefore necessary to use simulation tools which facilitate the formation and breaking of clusters, especially in isothermal−isobaric ensemble simulations that rely on translation and volume change moves to sample configuration space. In this work, we use AVBMC moves proposed by Chen and Siepmann to address these challenges.16 Two versions of AVBMC displacements are employed in our simulations. For single ions we use the AVBMC2 move, where a randomly selected ion is swapped between the bonded and nonbonded regions of another ion. The acceptance criteria used here are similar to those provided by Chen and Siepmann.16 These moves prove useful when sampling configurations consisting of charged clusters, which are known to appear in the vapor phase at certain conditions.4,15 To facilitate the sampling of configurations consisting of neutral clusters of various sizes we apply the AVBMC3 technique to displace ion-pairs. This algorithm allows growth of one cluster

(17)

k1 Nout

where is the number of ions of type i in the ‘out’ region of k1 (excluding ion k2) and Vin and Vout are the volumes of the ‘in’ and ‘out’ regions respectively. Similarly, the move (4) is accepted with the probability pacc

counter ⎫ ⎧ ⎪ Vout pj , n ⎪ πn Nink1 ⎨ ⎬ = min 1, counter ⎪ k1 ⎪ ⎩ πo (Nout + 1) Vin pj , o ⎭

(18)

We will show in Section 4 that these moves prove sufficient for sampling clusters that form within vapors consisting of ionic liquids. Phase Coexistence. To calculate liquid−vapor coexistence properties over a wide range of temperatures, we follow an approach outlined in a recent report from our group.2 Direct grand canonical simulations are first performed at select high temperatures. These direct simulations provide the density probability distribution over a range of densities that spans both the vapor and liquid regions, thereby providing the data necessary to locate a saturation point with the assistance of histogram reweighting.29 Saturation pressures psat (or activities ξsat) evaluated at these temperatures are used to formulate a guess of how these quantities evolve with temperature along the liquid−vapor saturation curve. In this work we assume a linear relationship between ln psat (or ln ξsat) and β to generate an initial guess. Isothermal−isobaric temperature expanded ensemble (NPT-TEE) or grand canonical temperature expanded ensemble (GC-TEE) simulations are then used to D



Article

One moves are partitioned such that the computational time invested in each move is roughly equal. The reservoir for each ion-type includes Nres = 500 ions in the ideal gas state. Before starting the simulation, we use a single molecular configuration to populate the reservoir. The reservoirs are initialized by completing 2000 to 5000 MC cycles, with a cycle defined as Nres MC steps. During a GC simulation, we generally complete two reservoir MC moves for every one MC move performed within the main simulation cell. The computational cost associated with maintaining the reservoir typically does not exceed 10% of the total computational effort. A partially grown (tagged) ion interacts with the rest of the system via a modified Lennard-Jones potential introduced by Lo and Palmer35

refine this guess over the course of a few iterative steps. Within the TEE simulations, subensembles are separated by a constant shift in β and the pressure (or activity) within each subensemble is set to match the current estimate for the saturation curve. NPT-TEE (or GC-TEE) subensemble probability distributions provide the relative Gibbs free energy Gk (or the grand potential Ωk) of a given phase k along the psat − β (or ξsat − β) path sampled. The psat − β (or ξsat − β) relationship is considered converged when the Gibbs free energy (or grand potential) of the liquid and vapor phases are sufficiently close at all temperatures of interest. Visited-states density probability distributions collected within each subensemble during TEE simulations and histogram reweighting techniques29 provide a means to deduce how Gk varies with p (or Ωk varies with ξ) at constant β, and thereby facilitate refinement of the saturation curve. For the ionic systems studied in this work, two to three iterations are required to converge to a solution. We typically perform relatively short runs at the outset of this process, and then complete longer runs with the converged pressures or activities. We note that the Gibbs free energy of the liquid phase changes by a statistically insignificant amount with the variation in psat during the iterative process, and therefore the NPT-TEE simulations of the liquid phase could be performed just once during this process.

⎧ ⎤12 ⎡ ⎤6 ⎫ ⎪⎡ ⎪ σ σ ⎬ u(r ) = λ 4ε⎨ − ⎢ ⎥ ⎢ ⎥ ⎪ ⎪ ⎣ r + ζ(1 − λ) ⎦ ⎭ ⎩⎣ r + ζ(1 − λ) ⎦ (19)

where r is the site−site separation distance, ε and σ are the conventional Lennard-Jones energy and size parameters, respectively, ζ is an adjustable parameter we set to 0.5σ, and λ = m/M is a scaling parameter that controls the growth of a molecule. Coulombic interactions are evaluated with the potential

3. FORCE FIELD AND SIMULATION DETAILS Molecular Models. We perform calculations with two ionic liquid models. The simpler model of the two is the restricted primitive model (RPM), where cations and anions are monovalently charged hard spheres of equal diameters. In what follows, the properties calculated for this model are scaled using hard sphere diameter σ as the size parameter and Q* = q2/(4πεεoσ) as the characteristic energy scale. Here, ε is the dielectric constant, εo is the permittivity of vacuum, and q is the absolute charge on the ions. For the realistic model, we use the classical united-atom force field for imidazolium-based ionic liquids introduced by Zhong et al.30 The force field was developed in the framework of AMBER/OPLS.31,32 The nonbonded interactions between pseudoatoms are described by Lennard-Jones 12−6 and point charges that interact via the Coulombic potential. It employs a harmonic potential for bond stretching and bond bending and a triple cosine potential to govern the sampling of dihedral and improper dihedral angles. This force field was selected because of its transferable nature which enables us to model imidazolium cations of varying alky chain lengths and different anions. Previous studies have reported that the simulations using this force field produce enthalpies of vaporization and critical points which are consistent with those deduced from experiments.3,4,30 We also complete a series of benchmarking calculations with propanol and benzene, for which we use the united-atom33 and explicit-hydrogen34 versions of the TraPPE model, respectively. Direct GC Simulations. A cubic simulation cell of volume 5 × 104 Å3 is used for [C1mim][BF4] and 1 × 105 Å3 is used for propanol and benzene. For the RPM, we use a cubic box of edge length L = 12. Additional calculations are performed with other cell sizes for [C1mim][BF4] and the RPM to examine finite-size effects. The MC move mix includes ion displacements and rotations, ion-pair displacements, HMC, HMC-One, as well as molecular insertions (or growth) and deletions (or reductions). Insertion/deletion moves are selected with 60% probability and the displacement, rotation, HMC, and HMC-

u(r ) =

λqiqj 1 4πε0 [r + ζ(1 − λ)]

(20)

where qi and qj represent the charges at sites i and j. Unlike [C1mim][BF4], insertions and deletions for the RPM are performed in a single step. A spherical cutoff coupled with a mean-field correction is used to account for long-range Lennard-Jones interactions.36 The Ewald summation36 with tinfoil boundary conditions is used to compute electrostatic interactions. Additional calculations are performed with vacuum boundary conditions for [C1mim][BF4] and the RPM to examine the influence of boundary conditions. Lennard-Jones interactions and the real space contribution to the Ewald summation are cutoff at a distance of 14 Å for [C1mim][BF4] and 16 Å for propanol and benzene. The algorithm described by Frenkel and Smit36 with a precision of 5 × 10−5 is used to determine the damping factor and number of lattice vectors associated with the Ewald summation. To reduce computational expense associated with the growth process we include the complete reciprocal and self-interaction contributions to the Ewald summation in the first growth stage and subsequently update the real space contribution only with a change in the stage number. GC-TEE Simulations. GC-TEE simulations are performed for the RPM and [C1mim][BF4] over a broad range of β. For the RPM, we complete calculations with simulation cells with edge length between L = 12 and 22 to facilitate the examination of system size effects. For [C1mim][BF4], we use a cubic simulation box of volume 1 × 105 Å over a range of β from 6.02 × 1019 J−1 (T ≈ 1203 K) to 1.208 × 1020 J−1 (T ≈ 600 K) with a discretization of Δβ = 2.0 × 1017 J−1. The reservoir temperature is updated to match the system temperature upon completion of a successful expanded ensemble move. As noted in our earlier work,2 the temperature dependence of zres is i determined via canonical temperature expanded ensemble simulations consisting of 100 isolated ions. We note that the E



Article

minimum temperature considered here is high relative to that used in NPT-TEE calculations due to difficulties encountered in insertion/deletion of complex molecules in the liquid phase at lower temperatures. Direct NPT Simulations. Direct NPT simulations are performed at a single temperature of T = 1100 K for [C1mim][BF4] with N = 100 molecules. Similar simulations at lower temperatures suffer from poor statistics due to difficulties associated with sampling inhomogeneous configurations that form at intermediate density. Direct NPT simulations are not effective for the RPM. The close association of ions within the vapor phase renders volume change moves involving the hard-sphere cores effective for very small volume changes only. This issue leads to slow sampling of density space. NPT-TEE Simulations. NPT-TEE simulations are performed for [C1mim][BF4] over a range of β from 6.02 × 1019 J−1 (T ≈ 1203K) to 2.654 × 1020 J−1 (T ≈ 273K) with a discretization of Δβ = 2.0 × 1017 J−1. A cubic simulation cell consisting of 100 ion pairs is employed for these simulations. Additional calculations are performed with 150 and 200 ion pairs to examine the effect of system size. The MC move mix includes ion displacements and rotations, ion-pair displacements, ion AVBMC2 displacements, ion-pair AVBMC3 displacements, HMC, HMC-One, and torsion-free HMC-One moves, volume change attempts and subensemble change moves. Volume changes are selected with a probability of approximately 1/N, subensemble change moves are attempted with probability 0.01, and the remaining moves are distributed using the “equal computational cost” approach noted above. The spherical cutoff distance for the Lennard-Jones and realspace Ewald contribution is taken to be a constant fraction of the simulation cell length, usually 50%. Computing Probability Distributions. To evaluate probability distributions (e.g., density or temperature probability distributions) we use a transition matrix Monte Carlo (TMMC) based scheme. The details of the overall strategy are provided in the appendix of a recent report from our group.2

Figure 1. GC density probability distributions for the RPM at T = 0.04294. The probability distribution Π(ρ) is normalized such that the integral over the entire density range is one. The density is defined as ρ = Nionsσ3/V. Blue and red lines indicate distributions obtained with tinfoil and vacuum boundary conditions, respectively. Solid, dashed, and dot−dashed lines correspond to the results obtained with cubic simulation boxes with L = 12, 18, and 22, respectively.

4. RESULTS AND DISCUSSION Figure 1 contains density probability distributions Π(ρ) for the RPM obtained from direct GC simulations at T = 0.04294 and liquid−vapor saturation conditions. We have plotted the logarithm of normalized probability distributions for system sizes of L = 12, 18, and 22 for both tinfoil and vacuum boundary conditions. Within these plots, the two peaks correspond to the coexisting vapor and liquid phases. We start by discussing the system size dependence of these profiles. It can be observed that the barrier between liquid and vapor phases increases with system size. This is due to the free energy costs involved in the generation of additional interfacial area. The probability distribution at intermediate densities can be used to extract important information about various interfacial configurations and to calculate surface tension.37−39 The figure also highlights the effect of boundary conditions used in the Ewald summation on liquid−vapor saturation properties. The profiles calculated by means of tinfoil and vacuum boundary conditions are statistically different for intermediate densities. This shows that the properties related to the interfacial configurations will be significantly affected by the choice of the boundary conditions used in the Ewald summation. Figure 2 shows the effect of system size and the nature of boundary conditions on liquid−vapor saturation properties

Figure 2. System size dependence of liquid−vapor coexistence properties for the RPM at T = 0.04294. Blue and red symbols correspond to properties obtained with tinfoil and vacuum boundary conditions, respectively. All quantities are plotted as a function of inverse volume. The quantity ln ξ∞ sat = −30.9856 denotes ln ξsat extrapolated to the macroscopic limit. The dashed lines simply connect adjacent data points.

computed from direct GC simulations at T = 0.04294. The data provided cover a range of box sizes spanning from L = 12 to 28. F



Article

We noted above the GC-TEE approach that we employ to obtain phase coexistence properties over a wide range of temperatures. Here, the temperature range is discretized into subensembles and the activity within each subensemble is refined over a series of iterative steps to deduce the ξsat(β) curve. Figure 3 provides data that illustrate how this process

At temperatures sufficiently below the critical temperature, bulk saturation properties are generally expected to scale linearly with reciprocal volume.40 We provide vapor pressures psat, liq vapor densities ρvap sat , liquid densities ρsat, saturated activities ln ξsat, and apparent surface tensions γlv as a function of reciprocal volume. We address the surface tension separately below, and now focus on the other four properties. Overall, we find that the properties calculated using vacuum boundary conditions show stronger system size dependence than those obtained via tinfoil boundary conditions. For ln ξsat and ρliq sat, both the tinfoil and vacuum data sets appear to scale linearly with V−1 and converge in the macroscopic limit. The tinfoil data are nearly independent of system size. In contrast, the vacuum data vary appreciably, with a 2% variation in ρliq sat and a shift of 0.02 in ln ξsat over the range of system sizes studied. For psat and ρvap sat , the tinfoil data sets appear to scale linearly with V−1, while, in contrast, the vacuum data sets show substantial curvature with V−1. Such curvature complicates extrapolation of the vacuum data to the infinite-system limit. As a result, it is difficult to appreciate how the tinfoil and vacuum results converge in the macroscopic limit. From a quantitative perspective, the tinfoilbased estimates for psat and ρvap sat vary by approximately 5% over the range of system sizes studied. Finally, we note that a similar system-size analysis was completed at T = 0.04790, yielding qualitatively similar results. Panagiotopoulos studied the dependence of system size on critical properties using GC simulation with the two boundary conditions discussed here.20 He showed that tinfoil- and vacuum-based critical properties converge in the infinite system size limit. Similar to the trends exhibited in Figure 2, he found that vacuum-based results showed stronger system size dependence than the analogous tinfoil results. Although the data for psat and ρvap sat raise some concerns, it generally appears that properties evaluated with tinfoil and vacuum boundary conditions converge in the macroscopic limit. Given the weaker system size dependence, we find tinfoil to be preferable to vacuum boundary conditions. This conclusion is consistent with the recommendations regarding this issue provided in earlier reports focused on ionic systems.36,41 We now return to the surface tension. The difference ΔΠ between the maximum (liquid, vapor peaks) and minimum (intermediate density) within a GC density probability distribution Π(ρ) (see Figure 1) is related to the finite-size surface tension γlv,L = kBTΔΠ/2L2.39 The finite-size scaling framework of Binder42 provides a means to extrapolate γlv,L data to the macroscopic limit. Specifically, he suggests that γlv,L scales with L−2. Given that the necessary data are available, we now take this opportunity to compute the surface tension using Binder’s approach. Using data gathered for the three largest system sizes (L = 18, 22, and 28), we obtain values of γlv = 0.00214(3) and 0.00213(2) with tinfoil and vacuum boundary conditions, respectively. Note that although γlv,L data are plotted versus L−3 in Figure 2, L−2 is used to complete the surface tension finite-size scaling analysis. For a point of comparison, we consider the work of Gonzalez-Melchor et al.,43 who determined the surface tension by analyzing the variation in the pressure tensor across an explicit liquid−vapor interface. These authors report a value of γlv = 0.0030(4) at T = 0.043, which is statistically different from the value we obtain using Binder’s finite-size scaling technique. Further investigation is necessary to understand the discrepancy between our value and that of Gonzalez-Melchor et al.

Figure 3. Convergence of RPM saturated activities with consecutive iterations of the GC-TEE method. Dashed lines indicate the fraction of the final saturated activity estimate ln ξfinal sat . Blue crosses correspond to the statistical uncertainty of ln ξfinal sat . Open red and green circles represent the absolute difference between ln ξfinal sat and ln ξsat obtained after the first and second iterations, respectively.

converges in just a few steps for the RPM. As a relevant point of reference, we plot the uncertainty in the saturated activity σln ξsat as a function of reciprocal temperature. These estimates are obtained from four independent simulation runs and provide a natural convergence limit. The dashed lines are guides that indicate levels of uncertainty relative to ln ξsat. As one might expect, σln ξsat evaluates to higher values at lower temperatures due to the difficulty in completing molecule transfers to and from the liquid phase at these conditions. Now consider Δln ξksat = |ln ξksat − ln ξfinal sat |, which is the difference between the k estimate of the saturated activity ln ξsat obtained after completion of iteration stage k and the final (converged) final 1 saturated activity ln ξsat . For example, the values ln ξsat correspond to those values computed using data generated with the linear relationship between ln ξ and β. We observe that Δln ξ1sat is in good agreement with σln ξsat at relatively high temperature, but is significantly larger than σln ξsat at moderate to low temperatures. However, after completion of the second iterative step, we find that the values for Δln ξ2sat track σln ξsat over the entire temperature range, indicating that the procedure converges within two iterations. Figures 4 and 5 provide the temperature dependence of saturated densities and vapor pressures, respectively, for the RPM. We provide results obtained with tinfoil boundary conditions from direct GC simulation (L = 12) and GC-TEE simulation (L = 12 and 22). Analogous direct GC simulation (L = 12) and GC-TEE simulation (L = 12) results obtained with vacuum boundary conditions are also provided. The figures also G



Article

most direct means for comparison, as these authors employed the same ensemble, system size, and boundary conditions as those used in one of the cases examined here, and we find excellent agreement. We also find reasonable agreement with the Gibbs ensemble results from Rai and Maginn.4 Our GC and GC-TEE results with similar system size and boundary condition are also consistent. Note that while this agreement is required at the single temperature used as a reference for GCTEE calculations (T = 0.04294 in our study), it is not guaranteed at other temperatures. The saturated vapor densities and vapor pressures obtained using the two different boundary conditions differ by a significant amount over the entire temperature range. In contrast, the corresponding saturated liquid densities show increasing compliance at lower temperatures. These results highlight the sensitivity of the vapor phase properties to the nature of the boundary condition. The GC-TEE results presented for the RPM liquid−vapor saturation curve extend to just above the expected triple point. At sufficiently low temperature the liquid freezes into a cesium chloride (CsCl) crystalline structure.44−46 Vega et al. estimate the vapor−liquid−solid (CsCl) triple point to be located at T = 0.0225.46 They arrive at this estimate by employing Gibbs− Duhem integration47 to construct the solid−liquid saturation curve and zero-pressure NPT simulations to estimate the density of the saturated liquid along the liquid−vapor saturation curve. Our results appear to be consistent with this estimate of the triple point, in the sense that we do not observe any signatures of crystal formation within our liquid-phase simulations. We now shift our attention to the united-atom model for [C1mim][BF4]. Figure 6 provides density probability distributions for [C1mim][BF4] collected at T = 1100 K from direct GC and NPT simulations. We again observe a substantial difference between distributions generated with vacuum and

Figure 4. Temperature dependence of the saturated densites for the RPM. Data are provided from direct GC simulations with tinfoil boundary conditions (blue circles) and vacuum boundary conditions (violet triangles) with L = 12. Solid black and dashed red lines correspond to data generated from GC-TEE simulations with tinfoil boundary conditions and L = 12 and 22, respectively. Dashed violet lines represent the data generated from GC-TEE simulations using vacuum boundary conditions and L = 12. Open green squares correspond to the grand canonical simulation data of Orkoulas and Panagiatopoulos24 employing vacuum boundary conditions and L = 12. Open cyan diamonds represent Gibbs ensemble data of Rai and Maginn4 with tinfoil boundary conditions. The black diamond represents the critical point reported by Orkoulas and Panagiatopoulos.24 The inset shows the liquid branch of the coexistence curve at the highest temperatures studied. Uncertainties for the TEE approaches are provided at select temperatures.

Figure 5. Temperature dependence of the vapor pressure for the RPM. Symbols are defined in the same manner as for Figure 4. The inset shows the portion of the plot at the highest temperatures studied.

Figure 6. GC (top) and NPT (bottom) density probability distributions for [C1mim][BF4] at T = 1100 K. The probability distribution Π(ρ) is normalized such that the integral over the entire density range is one. In the top figure the solid, dashed, dot-dashed, and dotted lines correspond to data obtained using a simulation box of volume 0.5 × 105 Å3, 1.0 × 105 Å3, 1.5 × 105 Å3, and 2.0 × 105 Å3, respectively. Blue and red lines indicate probability distributions obtained with tinfoil and vacuum boundary conditions, respectively.

contain Gibbs ensemble results from Rai and Maginn4 obtained with tinfoil boundary conditions and GC simulation (L = 12) results from Orkoulas and Panagiotopoulos24 obtained with vacuum boundary conditions. The latter study provides the H



Article

formation and breakup of these clusters. To probe the efficacy of these moves, we consider a series of simulations that are initiated with nonequilibrium configurations and follow a trajectory requiring the formation or breakup of clusters. We work at T = 600 K, where the equilibrated vapor consists of isolated ion pairs and a smaller proportion of neutral aggregates containing four and six ions. In one case, we start with a configuration containing relatively large aggregates and allow the system to relax using an isothermal−isobaric simulation. In the second case, we start with a configuration containing randomly distributed isolated ion pairs. Figure 8 contains the

tinfoil boundary conditions. However, these differences are small relative to those for the RPM. This contrast is reasonable given that the RPM is more ionic in nature than [C1mim][BF4]. The profiles obtained from NPT simulations exhibit more noise than those from GC simulations, especially at intermediate densities. These results point to the limitations of using direct NPT simulation to sample the density regime dominated by inhomogeneous fluid structures, an issue encountered with both ionic and nonionic systems.2 Due to the higher levels of uncertainty associated with direct NPT calculations, we use direct GC data to serve as a reference for both GC-TEE and NPT-TEE calculations. Figure 7 shows the system size dependence of vapor pressure, saturated densities, and saturated activities of

Figure 8. Evolution of configurational energy with the number of MC simulation cycles at T = 600 K. Blue and red curves correspond to the configurational energy from the NPT simulations initiated with isolated ion-pairs and large aggregates, respectively. The insets a and b show the initial and equilibrated cluster size distributions, respectively, for the simulation initiated with isolated ion-pairs. The insets c and d show the initial and equilibrated cluster size distributions, respectively, for the simulation initiated with large aggregates. The cluster size distribution is expressed as the fraction of the total number of aggregates.

Figure 7. System size dependence of liquid−vapor coexistence properties for [C1mim][BF4] at T = 1100 K. Symbols are defined in the same manner as in Figure 2.

[C1mim][BF4] calculated with direct GC simulation at T = 1100 K. Results for system sizes of 0.5 × 105 Å3, 1 × 105 Å3, 1.5 × 105 Å3, and 2 × 105 Å3 with tinfoil and vacuum boundary conditions are reported. While high uncertainties relative to the RPM case complicate the analysis, we find a number of similarities with the RPM. Specifically, the nature of the boundary condition has a clear, statistically significant impact on saturation properties, the tinfoil-based vapor pressures and densities are larger than the analogous vacuum-based estimates, and the tinfoil-based activity is statistically independent of system size while the vacuum-based activity is clearly systemsize dependent. From a quantitative perspective, the tinfoilbased estimates for the vapor pressure vary by approximately 10% over the range of system sizes studied, whereas the saturated liquid density is statistically independent of system size. A number of studies have highlighted the prominence of clusters within the vapor phase of ionic systems. These studies have shown that the proportion of larger aggregates increases with temperature and that isolated ion pairs are predominantly observed at lower temperatures.3,4 As noted above, we rely upon aggregation-volume-biasing moves16 to facilitate the

evolution of the configurational energy with the number of MC cycles for these two cases. The two trajectories converge to a common energy and density (not shown) following a few hundred MC cycles. The insets show the initial and final cluster size distributions calculated as the percentage of the total number of aggregates for the ‘breaking’ and ‘formation’ tests. We see that the final cluster size distributions for the two test cases are statistically equivalent. These tests provide confidence that the simulation tools employed in this work are able to sample the cluster size distribution of ionic liquids efficiently. In Section 2, we described the distance-biasing scheme we use for ion-pair transfers and displacements. To validate the approach used in our simulations, we compare the properties calculated from direct GC simulations of [C1mim][BF4] with and without distance-biasing moves. This study is performed at T = 1100 K (Tc < 1300 K) to allow acceptable sampling of states without distance-biasing. Validation tests are also performed with the non-ionic molecules propanol and benzene at T = 510 and 520 K, respectively. In these calculations, two molecules are transferred to or from the reservoir as a pseudo ion-pair. Table 1 provides the vapor pressures and saturated I



Article

Table 1. Saturation Properties Calculated with and without Using Distance-Biasing Moves Molecule

T (K)

psat (MPa) With bias

psat (MPa) Without bias

3 ρvap sat (kg/m ) With bias

3 ρvap sat (kg/m ) Without bias

3 ρliq sat (kg/m ) With bias

3 ρliq sat (kg/m ) Without bias

[C1mim][BF4] Propanol Benzene

1100 510 520

0.096(1) 3.397(7) 2.86(1)

0.096(1) 3.388(7) 2.863(3)

2.62(2) 89.2(4) 84.3(5)

2.64(2) 88.3(4) 84.4(1)

654(2) 491(3) 560(4)

655(2) 492(4) 562(3)

ature. The second window provides data over the relatively low temperature range of T = 273 to 333 K. Here, the values for Δp1sat are significantly larger than σpsat. After completion of the second iterative step, we find that Δp2sat is in better agreement with σpsat except for a few temperatures where p2sat deviates from final psat by more than 50%. Therefore, a third iteration is performed, with the resulting values for Δp3sat tracking σpsat over the entire temperature range, indicating that the process has converged. Figures 10 and 11 provide saturated densities and vapor pressures for [C1mim][BF4] over a range of temperatures

densities computed with and without the distance-biasing scheme. We find excellent agreement in all cases; ‘switching on/off’ distance-biasing does not impact saturation properties. The NPT-TEE approach provides an effective means for computing saturation properties over a wide range of temperatures. Figure 9 shows the convergence of vapor

Figure 9. Convergence of the [C1mim][BF4] vapor pressure at relatively high (top) and low (bottom) temperatures with consecutive iterations of the NPT-TEE scheme. Dashed black lines indicate percentages of the final vapor pressure estimate pfinal sat . Blue crosses correspond to the statistical uncertainty in pfinal sat . Open red, green, and violet circles represent the absolute difference between pfinal sat and the estimate for psat obtained after the first, second, and third iterations, respectively.

Figure 10. Temperature dependence of saturated densities for [C1mim][BF4]. Data are provided from direct GC (open blue circles), direct NPT (red triangles), GC-TEE (dashed blue lines), and NPTTEE (red lines) simulations. Open green squares correspond to the Gibbs ensemble data of Rai and Maginn.3 The inset shows the liquid branch of the coexistence curve at the highest temperatures studied. Uncertainties for the TEE approaches are provided at select temperatures.

pressures calculated from successive iterations of the NPTTEE technique. In a manner analogous to the GC-TEE analysis discussed above for the RPM, we plot the difference Δpksat = |pksat k − pfinal sat | between the estimate of the vapor pressure psat obtained after completion of iteration stage k and the final vapor pressure estimate pfinal sat along with the uncertainty in the vapor pressure σpsat. The first iteration is completed with a linear guess for the relationship between ln p and β, which is produced from direct GC data generated at temperatures of T = 1000 and 1100 K. The vapor pressure estimates stemming from this run are denoted as p1sat. These pressures are used to perform a second iteration, which results in p2sat, and so on. Due to the rather broad range of pressures traversed along the saturation line (almost 20 orders of magnitude), we analyze separately two temperature windows (note that just one NPT-TEE calculation is completed spanning the entire temperature range). The first window provides data over the relatively high temperature range of T = 667 to 1200 K, wherein we observe that values for Δp1sat are commensurate with σpsat, indicating that just a single iteration is sufficient to achieve convergence at high temper-

spanning from T = 273 to 1200 K. We present results obtained from direct GC, direct NPT, GC-TEE, and NPT-TEE simulations. We do not employ the iterative procedure described above to evaluate saturated activities for GC-TEE calculations. Instead, the activities are deduced from the converged vapor pressure estimates obtained from NPT-TEE simulations as explained in Appendix B. The GC-TEE simulations are terminated at relatively high temperature due to the difficulty of transferring ion-pairs to and from the dense liquid phase. Overall, we find very good agreement between the various methods employed here to compute saturation properties. We also show in Figures 10 and 11 the Gibbs ensemble results obtained by Rai and Maginn.3 The saturated liquid densities calculated from our simulations are in good agreement with their results. However, our saturated vapor J



Article

saturation properties at such low temperatures is to test the capabilities of the simulation strategies pursued in this work.

5. CONCLUSIONS In this paper we discussed a collection of Monte Carlo simulation strategies for calculating the liquid−vapor saturation properties of model ionic liquids. We outlined the use of distance-biasing and association-volume-biasing to address sampling difficulties encountered with ionic liquids due to their highly associating nature. Reservoir and growth expanded ensemble techniques were also discussed as a means to account for the intramolecular degrees of freedom of realistic models for ionic liquids. Several benchmarking calculations were completed with the RPM and a united-atom model for [C1mim][BF4] to examine the effects of system size and Ewald boundary condition. We also explored the efficacy of various techniques for computing saturation properties. We found that direct GC simulation provides an effective means to locate saturation points, particularly at relatively high temperatures along the coexistence curve. Direct NPT simulation proved less effective. The GC-TEE and NPT-TEE schemes provided a robust means to trace saturation curves to relatively low temperature. The GC-TEE approach was most effective for the RPM, as the hard nature of the potential coupled with strong ion association limited the efficacy of volume change moves. The NPT-TEE approach was most effective for the [C1mim][BF4] model, enabling us to compute saturation properties over a broad range of temperatures spanning from room temperature to the critical region. In general, we find the NPTbased approach provides the most effective means to trace saturation curves when working with force fields commonly used to describe realistic fluids.2 The system size and Ewald boundary condition employed to simulate ionic fluids have a significant impact on saturation properties. Considering system sizes typically employed in practice, these factors can account for differences of 20% or more in estimates for common properties, such as the vapor pressure. Overall, we found tinfoil boundary conditions provided a more robust means to compute thermodynamic properties. We note that this conclusion is consistent with earlier reports,20,36 wherein it has been argued that tinfoil boundary conditions are most appropriate for ionic systems. The strategies discussed in this work are general in nature and can be applied to study liquid−vapor phase behavior for a wide array of ionic fluids. For example, we are now applying these techniques to compute the saturation properties of model systems for which experimental data are available. The results generated from this effort will enable us to test the ability of common force fields to capture the saturation properties of RTILs. Finally, we note that the techniques pursued here are readily applicable to interfacial systems. In fact, we are now coupling the sampling strategies outlined here with interfacialpotential-based approaches1,49 to investigate the wetting properties of model ionic liquids.

Figure 11. Temperature dependence of the vapor pressure for [C1mim][BF4]. Symbols are defined in the same manner as in Figure 10. The violet diamonds correspond to the experimental data of Rocha et al.6 for [C2mim][NTf2]. The maroon cross indicates the vapor pressure psat = 1.70 × 10−17 MPa at T = 300 K. The inset shows the portion of the vapor curve at the highest temperatures studied.

densities and pressures are significantly lower than those from the Rai and Maginn study. At present, we cannot explain this discrepancy. We have explored a number of factors, but this process has not produced a satisfactory explanation. For example, the variations in vapor pressure due to ensemble, system size, and boundary condition, which are discussed above, cannot explain the discrepancy. We do note that the curvature of the Rai and Maginn vapor pressure data does not appear to be consistent with experimental48 and simulation (see Figure 5 and ref 2) data for a wide array of fluids. More specifically, in ln psat − 1/T space, the vapor pressure curve is generally concave, whereas the Rai and Maginn data form a convex curve. Before concluding we make a connection with the recent experimental study of Rocha et al.6 We include in Figure 11 their vapor pressure data for 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([C2mim][NTf2]) over the temperature range T = 445 to 483 K. Note that this is a different ionic liquid than we study here. Nonetheless, the comparison is interesting. The shape and slope of the experimental and simulation data are similar and the two curves are separated by just 1 order of magnitude, which is reasonable given the reported variation in vapor pressure with alkyl chain length and anion identity.3,4 The vapor pressure of ionic liquids at room temperature is of considerable interest. Our simulation data provide a value of 2 × 10−11 Pa for [C1mim][BF4] at 298 K. By way of comparison, the fitting equations provided by Rocha et al. give an extrapolated value of 4 × 10−10 Pa for [C2mim][NTf2] at 298 K. Again, the agreement is quite reasonable considering that the comparison involves different ionic liquids and that the vapor pressure data presented in Figure 11 span nearly 20 orders of magnitude. For completeness, we note that the melting point of pure [C1mim][BF4] at atmospheric conditions exceeds 373 K. Therefore, we are probing conditions that are not experimentally realizable. The motivation for calculating the

■

APPENDIX A: DETAILED ACCEPTANCE PROBABILITIES FOR GROWTH EXPANDED ENSEMBLE TRIAL MOVES We now return to the discussion related to the growth expanded ensemble and provide detailed working expressions for eqs 9 to 15 that we expect to be useful in practice. We consider moves that take the system from microstate o to K



Article

microstate n. We first consider the case in which complete ionpairs (M = 1) are transferred between the system and the reservoir. The acceptance probability for inserting a complete ion-pair consisting of ions of types i and j is given by

saturated activities at relatively low temperature from NPTTEE simulation data. This is of particular interest to our group, as saturated activities play an important role in methods we employ to compute interfacial properties.1 We now show how NPT-TEE simulation data are used to compute saturated activities ξsat. We take as input the saturated activity ξ1 at a single relatively high inverse temperature β1 (these are obtained via direct GC simulation) and the NPT-TEE subensemble probability distribution ΠNPT−TEE([β,p]) along the converged saturation curve. The saturated Gibbs free energy Gsat at a particular temperature is obtained from the probability distribution

⎧ ⎪ ⎛ V ⎞2 ⎟ pacc = min⎨1, ⎜ ⎪ ⎝ N + 1⎠ ⎩ res

ξe−β(En − Eo − Ei

⎫ ⎡ (N + 1)pcounter (rijcom) ⎤⎪ j,n ⎥⎬ ⎢⎣ W (rijcom)/W0(rijcom) ⎥⎦⎪ ⎭

− Ejres)⎢

(9′)

⎛ Π ([β , p]) ⎞ ⎟⎟ βGsat(β , p) = β1Gsat(β1 , p1 ) − ln⎜⎜ NPT − TEE ⎝ ΠNPT − TEE([β1 , p1 ]) ⎠

and the corresponding transition probability for removing a complete ion-pair is ⎧ ⎪ ⎛ N ⎞2 ⎛ 1 ⎞ pacc = min⎨1, ⎜ ⎟ ⎜ ⎟ ⎪ ⎝V ⎠ ⎝ξ⎠ ⎩ ⎡ W (r com)/W (r com) ⎤⎫ ⎪ ij 0 ij −β(En − Eo + Eires + Ejres)⎢ ⎥⎬ e ⎢⎣ Npcounter (rijcom) ⎥⎦⎪ j,o ⎭

+

i

where is the reservoir configurational partition function for an ion of type i at inverse temperature β. The ratio of reservoir configurational partition functions at temperatures β and β1 for a particular ion type (e.g., cation, anion) is obtained by performing a temperature expanded ensemble simulation with a collection of isolated ions that interact with the reservoir Hamiltonian.2 The saturated activity ξsat is computed as follows

(10′)

⎧ ⎡ ⎤⎫ ⎪ ⎪ res res 1 ⎥⎬ pacc = min⎨ 1, V 2ξe−β(En − Eo − Ei − Ej )⎢ com com ⎪ ⎪ ⎢⎣ W (rij )/W0(rij ) ⎥⎦⎭ ⎩

ln ξsat(β , p) = ln ξsat(β1 , p1 )

(11′)

and the corresponding transition probability for removing a tagged ion-pair is

+

⎧ ⎫ res res 1 = min⎨1, 2 e−β(En − Eo + Ei + Ej )[W (rijcom)/W0(rijcom)]⎬ ⎩ Vξ ⎭

■

(12′)

pacc = min{1, e

}

(13′)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: jerring@buffalo.edu. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge the financial support of the National Science Foundation (Grant No. CHE-1012356). Computational resources were provided in part by the University at Buffalo Center for Computational Research and the Rensselaer Polytechnic Institute Computational Center for Nanotechnology Innovations.

(14′)

and that of initiating destruction of a complete ion-pair (m = M in microstate o) is

■

ASSOCIATED CONTENT

The data presented in Figures 4, 5, 10, and 11 are provided. This material is available free of charge via the Internet at http://pubs.acs.org

⎫

⎧ ⎡ ⎤⎫ ⎪ ⎪ 1 pacc = min⎨1, N 2e−β(En − Eo)⎢ counter com ⎥⎬ ⎢⎣ Np ⎪ (rij ) ⎥⎦⎪ j,o ⎩ ⎭

(22)

* Supporting Information

⎧ ⎛ 1 ⎞2 −β(En − Eo) ⎟ e pacc = min⎨1, ⎜ [(N + 1) ⎝ ⎠ N 1 + ⎩ ⎭

N

S

The transition probability for completing the growth of a tagged ion-pair (m = M in microstate n) is

pjcounter (rijcom)]⎬ ,n

βGsat(β , p) − β1Gsat(β1 , p1 )

where N is the number of ion-pairs used to perform the NPTTEE simulation.

The acceptance probability for a stage change move wherein the molecule remains tagged is simply −β(En − Eo)

(21)

zres i (β)

The acceptance probability for inserting a tagged ion-pair is

pacc

⎛ z res(β) ⎞ i ⎟⎟ res ⎝ zi (β1) ⎠

∑ ln⎜⎜

■

(15′)

APPENDIX B: OBTAINING ACTIVITIES FROM NPT-TEE SIMULATION DATA In a previous report from our group,2 we showed that the NPTTEE approach provides advantages over the GC-TEE analog, particularly at low temperature where molecular exchanges become difficult to complete within the grand canonical scheme. It is therefore desirable to have a means to calculate

REFERENCES

(1) Rane, K. S.; Kumar, V.; Errington, J. R. J. Chem. Phys. 2011, 135, 234102−234113. (2) Rane, K. S.; Murali, S.; Errington, J. R. J. Chem. Theory Comput. 2013, 9, 2552−2566. (3) Rai, N.; Maginn, E. J. J. Phys. Chem. Lett. 2011, 2, 1439−1443. (4) Rai, N.; Maginn, E. J. Faraday Discuss. 2012, 154, 53−69. (5) Martin-Betancourt, M.; Romero-Enrique, J. M.; Rull, L. F. J. Phys. Chem. B 2009, 113, 9046−9049. L



Article

(41) Figueirido, F.; Del, B. G. S.; Levy, R. M. J. Chem. Phys. 1995, 103, 6133−6142. (42) Binder, K. Phys. Rev. A 1982, 25, 1699−1709. (43) Gonzalez-Melchor, M.; Alejandre, J.; Bresme, F. Phys. Rev. Lett. 2003, 90, 135506. (44) Smit, B.; Esselink, K.; Frenkel, D. Mol. Phys. 1996, 87, 159−166. (45) Vega, C.; Bresme, F.; Abascal, J. L. F. Phys. Rev. E 1996, 54, 2746−2760. (46) Vega, C.; Abascal, J. L. F.; McBride, C.; Bresme, F. J. Chem. Phys. 2003, 119, 964−971. (47) Kofke, D. A. Mol. Phys. 1993, 78, 1331−1336. (48) NIST Chemistry WebBook. http://webbook.nist.gov/chemistry/ (49) Grzelak, E. M.; Errington, J. R. J. Chem. Phys. 2008, 128, 014710.

(6) Rocha, M. A. A.; Lima, C. F. R. A. C.; Gomes, L. R.; Schroder, B.; Coutinho, J. A. P.; Marrucho, I. M.; Esperanca, J. M. S. S.; Rebelo, L. P. N.; Shimizu, K.; Lopes, J. N. C.; Santos, L. M. N. B. F. J. Phys. Chem. B 2011, 115, 10919−10926. (7) Zaitsau, D. H.; Fumino, K.; Emel’yanenko, V. N.; Yermalayeu, A. V.; Ludwig, R.; Verevkin, S. P. Chem. Phys. Chem. 2012, 13, 1868− 1876. (8) Zaitsau, D. H.; Kabo, G. J.; Strechan, A. A.; Paulechka, Y. U.; Tschersich, A.; Verevkin, S. P.; Heintz, A. J. Phys. Chem. A 2006, 110, 7303−7306. (9) Errington, J. R.; Panagiotopoulos, A. Z. J. Chem. Phys. 1999, 111, 9731−9738. (10) Duane, S.; Kennedy, A. D.; Pendleton, B. J.; Roweth, D. Phys. Lett. B 1987, 195, 216−222. (11) Escobedo, F. A.; de Pablo, J. J. J. Chem. Phys. 1996, 105, 4391− 4394. (12) Orkoulas, G.; Panagiotopoulos, A. Z. J. Chem. Phys. 1994, 101, 1452−1459. (13) Earle, M. J.; Esperanca, J. M. S. S.; Gilea, M. A.; Canongia Lopes, J. N.; Rebelo, L. P. N.; Magee, J. W.; Seddon, K. R.; Widegren, J. A. Nature 2006, 439, 831−834. (14) Leal, J. P.; Esperanca, J. M. S. S.; Da Piedade, M. E. M.; Lopes, J. N. C.; Rebelo, L. P. N.; Seddon, K. R. J. Phys. Chem. A 2007, 111, 6176−6182. (15) Chaban, V. V.; Prezhdo, O. V. J. Phys. Chem. Lett. 2012, 3, 1657−1662. (16) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 11275− 11282. (17) Binder, K. AIP Conf. Proc. 2003, 690, 74−84. (18) Challa, M. S. S.; Landau, D. P.; Binder, K. Phys. Rev. B 1986, 34, 1841−1852. (19) Wilding, N. B. Phys. Rev. E 1995, 52, 602. (20) Panagiotopoulos, A. Z. J. Chem. Phys. 2002, 116, 3007−3011. (21) Panagiotopoulos, A. Z. Fluid Phase Equilib. 1995, 104, 185−194. (22) Kristof, T.; Boda, D.; Liszi, J.; Henderson, D.; Carlson, E. Mol. Phys. 2003, 101, 1611−1616. (23) Yan, Q.; de Pablo, J. J. Phys. Rev. Lett. 2002, 88, 095504. (24) Orkoulas, G.; Panagiotopoulos, A. Z. J. Chem. Phys. 1999, 110, 1581−1590. (25) Hynninen, A.-P.; Panagiotopoulos, A. Z. Mol. Phys. 2008, 106, 2039−2051. (26) Caillol, J. M.; Levesque, D.; Weis, J. J. J. Chem. Phys. 2002, 116, 10794−10800. (27) Yan, Q.; de Pablo, J. J. J. Chem. Phys. 1999, 111, 9509−9516. (28) Hynninen, A.-P.; Dijkstra, M.; Panagiotopoulos, A. Z. J. Chem. Phys. 2005, 123, 084903. (29) Ferrenberg, A. M.; Swendsen, R. H. Phys. Rev. Lett. 1988, 61, 2635−2638. (30) Zhong, X.; Liu, Z.; Cao, D. J. Phys. Chem. B 2011, 115, 10027− 10040. (31) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765−784. (32) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. L. J. Comput. Chem. 1997, 18, 1955−1970. (33) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 3093−3104. (34) Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2007, 111, 10790− 10799. (35) Lo, C.; Palmer, B. J. Chem. Phys. 1995, 102, 925−931. (36) Frenkel, D.; Smit, B. Understanding Molecular Simulation; 2d ed.; Academic Press: London, 2001; Vol. 1. (37) MacDowell, L. G.; Shen, V. K.; Errington, J. R. J. Chem. Phys. 2006, 125, 034705−034715. (38) Schrader, M.; Virnau, P.; Binder, K. Phys. Rev. E 2009, 79, 061104. (39) Errington, J. R. Phys. Rev. E 2003, 67, 012102. (40) Landau, D. P.; Binder, K. A Guide to Monte Carlo Simulation in Statistical Physics; Cambridge University Press: Cambridge, 2000. M


Using Monte Carlo Simulation to Compute LiquidâVapor Saturation

Recommend Documents

Using Monte Carlo Simulation to Compute LiquidâVapor Saturation