Modeling of Mixing Acetone and Water - American Chemical Society

Apr 24, 2012 - Department of Inorganic and Analytical Chemistry, Budapest ... Group of Technical Analytical Chemistry, Szt. Gellért tér 4, H-1111 Bu...
0 downloads 0 Views 334KB Size
Article pubs.acs.org/JPCB

Modeling of Mixing Acetone and Water: How Can Their Full Miscibility Be Reproduced in Computer Simulations? Anita Pinke Department of Inorganic and Analytical Chemistry, Budapest University of Technology and Economics, Szt. Gellért tér 4, H-1111 Budapest, Hungary

Pál Jedlovszky* Laboratory of Interfaces and Nanosize Systems, Institute of Chemistry, Eötvös Loránd University, Pázmány P. Stny 1/A, H-1117 Budapest, Hungary, MTA-BME Research Group of Technical Analytical Chemistry, Szt. Gellért tér 4, H-1111 Budapest, Hungary, and EKF Department of Chemistry, Leányka u. 6, H-3300 Eger, Hungary ABSTRACT: The free energy of mixing of acetone and water is calculated at 298 K by means of thermodynamic integration considering combinations of three acetone and six water potentials. The Anisotropic United Atom 4 (AUA4) and Transferable Potential for phase Equilibria (TraPPE) models of acetone are found not to be miscible with any of the six water models considered, although the free energy cost of the mixing of any of these model pairs is very small, being below the mean kinetic energy of the molecules along one degree of freedom of 0.5RT. On the other hand, the combination of the Pereyra, Asar, and Carignano (PAC) acetone and TIP5P-E water models turns out to be indeed fully miscible, and it is able to reproduce the change of the energy, entropy, and Helmholtz free energy of mixing of the two neat components very accurately (i.e., within 0.8 kJ/mol, 2.5 J/(mol K), and 0.3 kJ/mol, respectively) in the entire composition range. The obtained results also suggest that the PAC model of acetone is likely to be fully miscible with other water models, at least with SPC and TIP4P, as well. properties of, among others, acetone−water mixtures,2,7,10,14,16 the problem of miscibility has only been realized two decades after the first acetone simulations. Weerasinghe and Smith realized that the reliable modeling of acetone−water mixtures requires either the simultaneous parametrization of the two models or the adjustment of the parameters of one model to that of the other ones. For this purpose, they modified the acetone model of the GROMOS96 force field26 by varying the fractional charges on the carbon and oxygen atoms of the carbonyl group to reproduce the experimental density of the equimolar mixture with SPC/E27 water.7 The obtained model, called the Kirkwood−Buff Force Field (KBFF), was then tested by comparing Kirkwood−Buff integrals as well as a number of other thermodynamic properties of its mixtures with SPC/E water with experimental data.7 In a subsequent study, Perera and Sokolić reported long molecular dynamics simulations of acetone−water mixtures, described by different combinations of acetone and water models.15 This study led to the striking conclusion that, with the exception of the KBFF−SPC/E

1. INTRODUCTION Acetone is one of the most widely used organic solvents, having numerous scientific and industrial applications. Most of these applications take advantage of the well-known fact that acetone is miscible with water in any proportion under ambient conditions. It is probably not as well-known, however, that the thermodynamic driving force behind this miscibility is very weak; the free energy change accompanying the mixing of the two neat components in any ratio does not exceed 0.65 kJ/mol at 298 K.1 This extremely small free energy gain of their mixing led to severe difficulties in modeling acetone−water mixtures by computer simulation in the past two decades. Due to its importance in several fields of science, acetone was one of the first molecular liquids studied also by computer simulation methods. Since the first simulation of liquid acetone, reported by Evans and Evans in 1983,2 a number of various potential models have been developed2−10 and used to simulate the properties of neat liquid acetone2,5,11,12 as well as that of aqueous acetone solutions,3,7,10,12−18 nonaqueous19,20 and ternary mixtures of acetone,21 the liquid−vapor interface of neat acetone,22,23 and the adsorption layer of acetone on ice.24,25 Because these models were, in general, able to satisfactorily reproduce a number of different experimental © 2012 American Chemical Society

Received: March 19, 2012 Revised: April 18, 2012 Published: April 24, 2012 5977

dx.doi.org/10.1021/jp302629r | J. Phys. Chem. B 2012, 116, 5977−5984

The Journal of Physical Chemistry B

Article

SPC/E,27 the four-site TIP4P30 and TIP4P/2005,34 and the five-site TIP5P35 and TIP5P-E36 ones in combination with both acetone potentials. Since, to our opinion, a definite answer to the question of miscibility cannot be given without accurate and extensive free energy calculations, here we determine the free energy change accompanying the mixing of the two models considered in the entire composition range using the method of thermodynamic integration.28,29 During the course of this work, Pereyra, Asar, and Carignano reported computer simulations of acetone−water mixtures using a reparametrized version of the CHARMM27 acetone model37,38 and TIP5P-E water.10 The main conceptual development behind the model of Pereyra, Asar, and Carignano (referred to here as the PAC model) is that the magnitude of the fractional charges depends on the composition of the system.10 Although this treatment restricts the use of the PAC model to acetone−water binary mixtures, its success in accurately reproducing the experimental heat of mixing in the entire composition range stresses the importance of accounting for the polarizability of the acetone molecule when reproducing the thermodynamic properties of acetone−water mixtures. To test also this concept we included this acetone model in combination with TIP5P-E water in the present investigation. The paper is organized as follows. In Section 2 the tested acetone and water models are described. The method of thermodynamic integration is briefly reminded in Section 3, while in Section 4 details of the calculations performed are given. The obtained results are presented and discussed in detail in Section 5. Finally, in Section 6 the main conclusions of this study are summarized.

model pair, all model combinations show phase separation at least at certain compositions. This finding provoked the question why this phase separation has never been noticed before. The reason for this is certainly related to the fact that the free energy difference between the two neat liquids and their mixture is very small even for nonmiscible model pairs, and hence demixing occurs only in longer time scales than what was used in bulk phase simulations in the 90s. Thus, for instance, Perera and Sokolić observed demixing of SPC/E water and OPLS acetone4 only after a 2 ns long equilibration, while after an equilibration period of 0.5 ns the system still looked homogeneous (see Figure 7 of ref 15). Observation of demixing is further hindered by the use of periodic boundary conditions and small system sizes, making distinguishing between phase separation and self-association difficult. Recently, we performed extensive free energy calculations using the method of thermodynamic integration28,29 to study the miscibility of the rigid version of the KBFF acetone model with both SPC/E and TIP4P30 water in the low acetone mole fraction range18 and reached the conclusion that the free energy changes accompanying the mixing of these model pairs are, although slightly, positive. This finding is in contradiction with the former results of both Weerasinghe and Smith (using the original, flexible KBFF model)7 and Perera and Sokolić,15 who both claimed miscibility of KBFF acetone with SPC/E water. The conclusion of the latter authors was based only on a simple visual inspection of snapshots of the system simulated, and hence, strictly speaking, they only concluded that during the course of their simulation they did not observe demixing. Weerasinghe and Smith claimed miscibility on the basis of calculating the Kirkwood−Buff integrals of the system, without performing calculation of the free energy of mixing itself. Unfortunately, Kirkwood−Buff integrals can be calculated in computer simulations with rather large uncertainties which can lead to dubious results when the effect is too small. The question of what is the most accurate way of detecting demixing is still under debate;31,32 the importance of this issue is emphasized by the fact that clearly the thermodynamic driving force of mixing (in the work of Weerasinghe and Smith7) or demixing (in our previous work18) of the two components is found to be very small. Furthermore, the possibility that the observed qualitatively different mixing behavior originates in the different treatment of the internal degrees of freedom (i.e., the use of rigid vs flexible acetone model) cannot be excluded either. Nevertheless, no acetone− water model pair has been undoubtedly shown to be fully miscible so far. In this paper, we aim at investigating the problem of miscibility using two relatively new acetone models. We have chosen the potentials belonging to the Transferable Potential for phase Equilibria (TraPPE)8 and Anisotropic United Atom 4 (AUA4)9 force fields as promising candidates in this respect. The TraPPE model is a reparametrized version of the OPLS force field, aiming at a better reproduction of various phase equilibria, while the basic idea behind the AUA4 model is that moving the center of the Lennard-Jones interaction from the C atom toward the H atoms of the CH3 group (treated as a united atom) improves the reproduction of a set of various, including thermodynamic, properties of the system. To avoid strong dependence of the obtained results on the particular choice of the water model and to increase the chance of finding a fully miscible acetone−water potential pair, we have used six different water models, namely, the three-site SPC33 and

2. POTENTIAL MODELS All the potential models used in this study are rigid and nonpolarizable, and the energy of their interaction is pairwise additive; i.e., the total energy of the system can be calculated as the sum of the interaction energies of the molecule pairs. The intermolecular interactions are described by Coulomb terms acting between the fractional charges and Lennard-Jones interactions. Hence, the uij interaction energy of molecules i and j can be calculated as ni

uij =

nj

∑∑ A=1 B=1

12 ⎡⎛ ⎛ σ ⎞6 ⎤ σAB ⎞ 1 qAqB ⎢ ⎟⎟ − ⎜⎜ AB ⎟⎟ ⎥ + 4εAB⎢⎜⎜ ⎥ 4πε0 riA , jB r ⎝ riA , jB ⎠ ⎦ ⎣⎝ iA , jB ⎠

(1)

Here the indices A and B run over the ni interaction sites of molecule i and nj interaction sites of molecule j, respectively; qA and qB are the fractional charges carried by the respective interaction sites; ε0 is the vacuum permittivity; σAB and εAB are the Lennard-Jones distance and energy parameters, respectively, calculated from the values corresponding to the individual sites using the Lorentz−Berthelot rule,39 i.e. σAB =

σA + σB 2

(2)

εAεB

(3)

and

εAB =

and riA,jB is the distance between site A of molecule i and site B of molecule j. 2.1. Acetone Models. In this study, we considered three acetone models. In the TraPPE model of acetone,8 the 5978

dx.doi.org/10.1021/jp302629r | J. Phys. Chem. B 2012, 116, 5977−5984

The Journal of Physical Chemistry B

Article

models. In all of these models, one single Lennard-Jones interaction center is placed on the O atom, and qH positive fractional charges are carried by the two H atoms. In the three-site SPC33 and SPC/E27 models, the compensating −2qH fractional negative charge is also located on the O atom, while in the four-site TIP4P30 and TIP4P/200534 models this charge is displaced from the O atom along the main symmetry axis of the molecule toward the H atoms to the nonatomic site M. In the five-site TIP5P35 and TIP5P-E36 models, the positive charges carried by the H atoms are compensated by two negative fractional charges of −qH, carried by the nonatomic interaction sites L, placed in tetrahedral directions around the O atom with respect to the two H atoms to mimic the lone electron pairs. The interaction and geometry parameters of these water models are summarized in Tables 3 and 4, respectively.

Lennard-Jones and Coulomb interactions are centered on the atomic sites, and the methyl groups are treated as united atoms. The AUA4 model9 also treats the CH3 groups as united atoms; however, the center of their Lennard-Jones interaction is displaced from the C atom by the distance δ toward the centerof-mass of the three H atoms (i.e., along the extension of the C−C bond) to account also for the nonspherical shape of the CH3 group. All the other interactions are centered on the respective atomic sites also in this model. The PAC model10 treats all the atoms of the acetone molecule as separate interaction centers and contains no nonatomic interaction sites. Here we used a rigid version of the originally flexible model, fixing the bond lengths and bond angles to their equilibrium values. This model, developed specifically for acetone−water mixtures, accounts for the polarization of the acetone molecules due to their strongly polar water neighbors in an average way, as all the fractional charges of the model depend on the mole fraction of acetone in its aqueous solution simulated, xac, as10 q(xac) = q(1)(1.1502 − 0.2385xac +

0.0883xac2)

Table 3. Interaction Parameters of the Water Models Used

(4)

where the fractional charges corresponding to xac = 1, i.e., to neat acetone, q(1), along with all other parameters of the model are taken from the CHARMM27 force field.37,38 The interaction parameters as well as bond lengths and bond angles of these models are collected in Tables 1 and 2, respectively.

site

σ/Å

ε/kJ mol−1

TraPPE

CH3 C O CH3 C O H C C(O) O

3.79 3.82 3.05 3.607 3.020 2.981 2.352 3.671 3.564 3.029

0.8144 0.3324 0.6565 0.9985 0.5143 0.8020 0.0920 0.3347 0.2929 0.5020

AUA4

PAC

q/e

σ/Å

ε/kJ mol−1

qH/e

SPC SPC/E TIP4P TIP4P/2005 TIP5P TIP5P-E

3.166 3.166 3.154 3.1589 3.120 3.097

0.6502 0.6502 0.6491 0.7749 0.6699 0.7453

0.41 0.4238 0.52 0.5564 0.241 0.241

Table 4. Geometry Parameters of the Water Models Used

Table 1. Interaction Parameters of the Acetone Models Used model

model

model

δ/Å

0.424 −0.424 0.216 0.49 −0.49 0.09a −0.27a 0.55a −0.55a

a

site pair

distance (Å)

SPC

O−H

1.000

SPC/E

O−H

1.000

TIP4P

O−H O−M

0.9572 0.15

TIP4P/2005

O−H O−M

0.9572 0.1546

TIP5P

O−H O−L

0.9572 0.70

Values refer to neat acetone; in aqueous mixtures the charges scale according to eq 4.

Table 2. Geometry Parameters of the Acetone Models Used model

bond

bond length (Å)

TraPPE

CO C−CH3

1.229 1.520

AUA4

PAC

CO C−CH3

CO C−C C−H

angle

bond angle (deg)

CH3−CO CH3−C−CH3

121.4 117.2

CH3−CO CH3−C−CH3

120.4 119.2

C−CO C−C−C H−C−C H−C−H

122.0 116.0 110.5 108.4

TIP5P-E

O−H O−L

angle formed by

angle (deg)

H−O−H

109.47

H−O−H

109.47

H−O−H

104.52

H−O−H

104.52

H−O−H L−O−L

104.52 109.47

H−O−H L−O−L

104.52 109.47

0.9572 0.70

3. CALCULATION OF THE FREE ENERGY OF MIXING AND THE METHOD OF THERMODYNAMIC INTEGRATION The method of thermodynamic integration28,29 provides the difference of the Helmholtz free energy, A, between the state of interest, denoted as Y, and an appropriately chosen reference state, X. The free energy difference is calculated as an integral over a continuous path connecting the state of interest and the reference state through a set of fictitious intermediate states

1.229 1.522

1.230 1.522 1.111

ΔA = AY − A x =

∫0

1

∂A(λ) dλ ∂λ

(5)

Here the coupling parameter λ defines the continuous path along which the system is moved from state X to Y. The value

2.2. Water Models. In this study, we considered two widely used representatives of the three-, four-, and five-site water 5979

dx.doi.org/10.1021/jp302629r | J. Phys. Chem. B 2012, 116, 5977−5984

The Journal of Physical Chemistry B

Article

Having the integrand of eq 13 evaluated at a number of λ values, the excess Helmholtz free energy can be obtained by performing the integration numerically. Given that a simulation is also performed in state Y, i.e., at λ = 1, where the internal energy of the system can also be calculated, the excess entropy can also be evaluated as

of λ changes between 0 and 1, being λ = 0 in the reference state, X, and λ = 1 in the state of interest, Y. Using the fundamental statistical mechanical equations A = −kBT ln Q

(6)

and Q=

∫ exp(−βU )dq N

S=

(7)

the integrand of eq 5 can be written as −kBT ∂Q ∂A = N ∂λ ∫ exp(−βU )dq ∂λ

In these equations, Q is the canonical partition function; U is the potential energy; T is the absolute temperature; kB is the Boltzmann constant; β = 1/kBT; and qN stands for the 3N position coordinates of the N particles; i.e., it represents the position of the system in the configurational space. Performing the derivation in eq 8 using also eq 7, one obtains ∂U

∂U (λ) ∂λ

λ

(9)

where the brackets ⟨....⟩λ denote ensemble averaging at a given λ value. The U(λ) function, i.e., the energy of the fictitious intermediate states along the continuous transition path, is usually defined by the energy of the state of interest, UY, and that of the reference state, UX, as U (λ) = λ k UY + (1 − λ)k UX

(10)

A mix = ΔA sol − xacΔA ac − (1 − xac)ΔA wat

where the exponent k is usually chosen to be 4 since in threedimensional systems in which the intermolecular repulsion decays with r−12 the use of smaller exponents leads to a singularity at the λ = 0 end of the integral of eq 5.28 Equation 10 satisfies the boundary conditions of U(0) ≡ UX and U(1) ≡ UY. When the excess Helmholtz free energy of a condensed phase relative to the ideal gas state of the system is calculated, the energy of the reference state is UX ≡ 0, and hence eq 10 simplifies to U (λ) = λ k UY

+ RT[xac ln xac + (1 − xac)ln(1 − xac)] Umix = ΔUsol − xacΔUac − (1 − xac)ΔUwat

− R[xac ln xac + (1 − xac)ln(1 − xac)]

4. COMPUTATIONAL DETAILS The integrand of eq 13 has been evaluated at six different λ values by performing Monte Carlo simulations. Five of the used λ values, being 0.046911, 0.230765, 0.5, 0.769235, and 0.953089, correspond to the five point Gaussian quadrature, while the sixth λ value has been chosen to be 1 to evaluate the energy and hence also the excess entropy of the systems. The integration has been performed analytically, by fitting a fifthorder polynomial to the integrand values obtained at the six λ points. The integrand values calculated at the six λ points considered together with the fitted polynomials are shown in Figure 1 for the systems containing TIP5P-E water at the acetone mole fraction values of 0.2, 0.5, 0.7, and 1.0 for the three acetone models considered here. To evaluate the λ ensemble average in eq 13, Monte Carlo simulation has been performed on the canonical (N,V,T) ensemble at the T* virtual temperature corresponding to the actual λ value. The real temperature of the system, T, has been

Substituting this expression to eq 5, the excess Helmholtz free energy of state Y relative to the ideal gas state of the system can be calculated as 1

λ k − 1 < UY>λ dλ

(18)

respectively. In these equations, the subscripts sol, ac, and wat refer to the acetone−water mixture, neat acetone, and neat water, respectively, and Δ refers to the excess of the corresponding quantity with respect to the ideal gas state.

(12)

∫0

(17)

Smix = ΔSsol − xacΔSac − (1 − xac)ΔSwat

and eq 9 can be rewritten as

ΔA = AY − A id.gas = k

(16)

and

(11)

∂A = kλ k − 1 λ ∂λ

(15)

To calculate the Helmholtz free energy of mixing of two neat components, we proposed to use the following thermodynamic cycle.40 The two neat systems are first brought to the ideal gas state; the Helmholtz free energy change accompanying this step can be calculated by thermodynamic integration. Then the two components are mixed in the ideal gas state. This step is accompanied by the Helmholtz free energy change of the ideal mixing, i.e., RT(x1 ln x1 + x2 ln x2), where x1 and x2 are the mole fractions of the two components. Finally, the mixed system is brought back from the ideal gas state to the state of interest; the corresponding Helmholtz free energy change is again accessible by thermodynamic integration (see Figure 1 of ref 40). (It should be noted that this thermodynamic cycle misses the term corresponding to the volume change of the neat components in the liquid state upon mixing; however, the free energy change associated with this step is negligibly small in the case of systems of liquid-like densities.) Considering this thermodynamic cycle and using eq 15, the full molar Helmholtz free energy, internal energy, and entropy change accompanying the mixing of neat acetone and water can be evaluated as

(8)

∫ exp(−βU )dq N ∂A = ∂λ = ∂λ ∫ exp(−βU )dq N

U−A T

(13)

In evaluating the integrand of eq 13, a computer simulation has to be performed on the canonical ensemble at the state λ (i.e., using the potential function of eq 11). However, the Boltzmann factor used in the ensemble averaging at λ can be rewritten as exp( −U (λ)/kBT ) = exp(−λ k UY /kBT ) = exp(−UY /kBT *) (14)

where T* = T/λk. Therefore, the simulation performed at the real temperature of the system, T, with the potential function U(λ) is technically equivalent with a simulation performed at the virtual temperature T* with the full potential function, UY. 5980

dx.doi.org/10.1021/jp302629r | J. Phys. Chem. B 2012, 116, 5977−5984

The Journal of Physical Chemistry B

Article

number of simulations performed in this study has been above 750. The simulations have been performed using the program MMC.42 In each simulation step, a randomly chosen molecule has been randomly translated and randomly rotated around a randomly chosen space-fixed axis by no more than 0.25 Å and 20°, respectively. At least 25% of the attempted moves have been successful in every case; at high virtual temperatures (i.e., close to the ideal gas state), the rate of accepted and tried moves has been close to 1. All intermolecular interactions have been truncated to zero beyond the center−center cutoff distance of 12.5 Å, and the long-range part of the Coulomb interaction has been accounted for by the method of reaction field correction39,43,44 under conducting boundary conditions. The systems have been equilibrated in a 5 × 107 Monte Carlo steps long run at every T* virtual temperature (i.e., at every λ value). After equilibrium is reached, the ensemble average in eq 13 has been evaluated over a 108 Monte Carlo steps long trajectory in every case.

5. RESULTS AND DISCUSSION The free energy of mixing of the AUA4 and TraPPE acetone models with the six water models considered is shown in Figure 2 in the entire composition range. The uncertainty of the obtained values does not exceed ±0.1 kJ/mol in any case. As is seen, all model combinations result in positive Amix values at every composition, with the exception of AUA4 acetone and TIP4P water at xac = 0.9. Although this finding clearly indicates that none of these model combinations can form thermodynamically stable mixtures with each other, the free energy cost of their mixing is very small, being below 0.5RT, i.e., the average kinetic energy of the molecules along one degree of freedom in almost every case. The obtained small values of Amix clearly stress the difficulty of correctly describing the mixing behavior of acetone with water even qualitatively. It is also seen that the AUA4 acetone model gives lower free energy of mixing values than TraPPE with any of the six water models considered. Further, the Amix values obtained with the SPC, TIP4P, and TIP5P-E water models are consistently lower than those with SPC/E, TIP4P/2005, and TIP5P, respectively, considering mixing with either of the two acetone models. These findings indicate that both the decrease of the water and the increase of the acetone dipole moment lead to smaller Amix values, i.e., bring it closer to zero. The free energy, internal energy, and entropy values accompanying the mixing of the PAC acetone model with TIP5P-E water are shown and compared with the experimental data1 in Figure 3. Clearly, unlike all other model combinations considered, this model pair results in negative Amix values in the entire composition range. It is also seen that the calculated data follow reasonably well the shape of the experimental curve for all three mixing quantities. Furthermore, all the calculated values agree excellently with the experimental data, and their deviation never exceeds 0.3 kJ/mol for Amix, 0.8 kJ/mol for Umix, and 2.5 J/(mol K) for Smix (i.e., 0.75 kJ/mol for TSmix); i.e., these differences remain well below the average energy of the thermal motion along one degree of freedom of 0.5RT. (When comparing the calculated values with the experimental data, it should be noted that the experimental data refer to constant-pressure conditions; i.e., they correspond to the Gibbs free energy and enthalpy rather than to the Helmholtz free energy and energy of mixing of the two components. However, these quantities differ only in the PV term, which is negligibly

Figure 1. Integrand of the thermodynamic integration for acetone− water mixtures with the acetone mole fractions of 0.2 (squares), 0.5 (circles), 0.7 (up triangles), and 1.0 (down triangles) calculated at the selected λ values in the systems containing TIP5P-E water and AUA4 (top panel), TraPPE (middle panel), and PAC (bottom panel) acetone. The polynomials fitted to these data are shown by full lines.

298 K. The cubic basic simulation box has contained 864 molecules; the length of its edges has been set to reproduce the experimental density of the system.41 Standard periodic boundary conditions have been applied. Mixtures of nine different compositions have been simulated; the acetone mole fraction of these mixed systems, xac, has been changed between 0.1 and 0.9 with an increment of 0.1. In addition, to evaluate the Helmholtz free energy of mixing according to eq 16, simulations of the two neat systems have also been performed. The number of the acetone and water molecules (Nac and Nwat, respectively) as well as the edge length L of the basic simulation box of the systems of different compositions simulated are collected in Table 5. Table 5. Characteristics of the Acetone−Water Systems Simulated xac

Nac

Nwat

L/Å

0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00

0 86 172 259 346 432 518 605 692 778 864

864 778 692 605 518 432 346 259 173 86 0

29.57 32.11 34.40 36.50 38.25 40.02 41.66 43.22 44.67 46.05 46.89

We have investigated all 12 possible combinations of the TraPPE and AUA4 acetone models with the six water models considered. In addition, following the idea of Pereyra et al.10 the mixture of the PAC model of acetone with TIP5P-E water has also been studied. Considering that the thermodynamic integration required six simulations per system, the total 5981

dx.doi.org/10.1021/jp302629r | J. Phys. Chem. B 2012, 116, 5977−5984

The Journal of Physical Chemistry B

Article

Figure 2. Free energy of mixing of the AUA4 (left) and TraPPE (right) models of acetone with the SPC (black squares), SPC/E (red circles), TIP4P (green up triangles), TIP4P/2005 (dark blue down triangles), TIP5P (light blue diamonds), and TIP5P-E (magenta side triangles) models of water as a function of the acetone mole fraction. The middle bar shows the mean kinetic energy of the molecules along one degree of freedom, 0.5RT.

components, as it only indicates stability relative to the two neat systems, but not to two mixtures of different compositions. Two components form a thermodynamically stable binary mixture on the canonical ensemble if and only if the following condition is satisfied:45 D = 1 + xac(1 − xac)

∂ 2(A mix /RT ) ∂ 2xac

>0 (19)

Therefore, to test the full miscibility of the TIP5P-E water and PAC acetone models, the value of D also has to be determined in the entire composition range. Here we used two different ways for evaluating D. First, we performed double direct derivation of the calculated Amix(xac) data. Second, we fitted the Margules equation1 (A mix /kBT ) = xac(1 − xac)[A12 xac + A 21(1 − xac) − (C12xac + C21(1 − xac))xac(1 − xac) + Bxac2(1 − xac)2 ]

(20)

to these data and performed the double derivation analytically for the fitted function. The fitting of eq 20 to the Amix(xac) data resulted in the parameter values of A12 = −1.513, A21 = −2.600, C12 = −4.174, C21 = −4.416, and B = −4.448; the fitted curve is also shown in Figure 3. The D values calculated in both ways are shown in the entire composition range in Figure 4. The data obtained by direct numerical double differentiation scatter around the D(xac) curve derived from the Margules fit of the data. As is seen, the D values obtained in both ways clearly exceed zero and hence fulfill the inequality of eq 19 at any acetone mole fraction value. This result provides clear evidence that, similarly to real acetone and water, the PAC acetone and TIP5P-E water models are indeed miscible with each other at any ratio.

Figure 3. Helmholtz free energy (top panel), internal energy (middle panel), and entropy (bottom panel) of mixing of the PAC acetone and TIP5P-E water models as a function of the acetone mole fraction (full circles). Experimental data (ref 1) are shown by solid lines; the Margules equation (see eq 20) fitted to the simulated Amix(xac) data is indicated by a dashed line.

small, being in the order of a few joules/mole for liquid systems at atmospheric pressure.) It is also seen that the calculated energy and entropy of mixing values deviate somewhat stronger from the experimental data than Amix, in particular, at higher acetone mole fraction values, indicating that, besides the excellent performance of the model pair in this respect, an energy−entropy compensation effect also helps in accurately reproducing the free energy of mixing. It should also be noted that the enthalpy of mixing values reported by Pereyra et al. for this model pair agree even considerably better with the experimental data than the energy of mixing values reported here; the 0.1−0.5 kJ/mol difference between the two sets of calculated values can be attributed to the flexible vs rigid treatment of the acetone molecules. A negative value of the free energy of mixing is a necessary but not sufficient condition of the miscibility of two

6. SUMMARY AND CONCLUSIONS In this paper, we presented extensive calculations of the free energy of mixing of acetone and water, considering three acetone and six water models, by means of the method of thermodynamic integration. The obtained results showed that, similarly to other acetone models tested previously,15,18 neither the AUA4 nor the TraPPE model is miscible with water as described by any of the six models considered here. Comparison of the free energies of mixing obtained with 5982

dx.doi.org/10.1021/jp302629r | J. Phys. Chem. B 2012, 116, 5977−5984

The Journal of Physical Chemistry B



Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This project is supported by the Hungarian OTKA Foundation under project No. 75328. The authors are indebted to Dr Gábor Jancsó (KFKI AEKI, Budapest, Hungary) for critically reading the manuscript.



REFERENCES

(1) Villamañań , M. A.; Van Ness, H. C. J. Chem. Eng. Data 1984, 29, 429−431. (2) Evans, G. J.; Evans, M. W. J. Chem. Soc., Faraday Trans. II 1983, 79, 153−165. (3) Ferrario, M.; Haughney, M.; McDonald, I. R.; Klein, M. L. J. Chem. Phys. 1990, 93, 5156−5166. (4) Jorgensen, W. L.; Briggs, J. M.; Contreras, M. L. J. Phys. Chem. 1990, 94, 1683−1686. (5) Jedlovszky, P.; Pálinkás, G. Mol. Phys. 1995, 84, 217−233. (6) Hermida-Ramón, J. M.; Ríos, M. A. J. Phys. Chem. A 1998, 102, 2594−2602. (7) Weerasinghe, S.; Smith, P. E. J. Chem. Phys. 2003, 118, 10663− 10670. (8) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2004, 108, 17596−17605. (9) Ferrando, N.; Lachet, V.; Boutin, A. J. Phys. Chem. B 2010, 114, 8680−8688. (10) Pereyra, R. G.; Asar, M. L.; Carignano, M. A. Chem. Phys. Lett. 2011, 507, 240−243. (11) Bushuev, Yu. G.; Davletbaeva, S. V. Russ. Chem. Bull. 1999, 48, 25−34. (12) Gomide Freitas, L. C.; Cordeiro, J. M. M.; Garbujo, F. L. L. J. Mol. Liq. 1999, 79, 1−15. (13) Bushuev, Yu. G.; Korolev, V. P. Russ. Chem. Bull. 1998, 47, 569− 577. (14) Idrissi, A.; Longelin, S.; Sokolić, F. J. Phys. Chem. B 2001, 105, 6004−6009. (15) Perera, A.; Sokolić, F. J. Chem. Phys. 2004, 121, 11272−11282. (16) Idrissi, A; Damay, P. J. Non-Cryst. Solids 2006, 352, 4486−4489. (17) Jedlovszky, P.; Idrissi, A. J. Chem. Phys. 2008, 129, 164501−1−7. (18) Jedlovszky, P.; Idrissi, A.; Jancsó, G. J. Chem. Phys. 2009, 130, 124516−1−8. (19) Venables, D. S.; Schmuttenmaer, A. C. J. Chem. Phys. 2000, 113, 3249−3260. (20) Idrissi, A.; Vyalov, I.; Kiselev, M.; Jedlovszky, P. Phys. Chem. Chem. Phys. 2011, 13, 16272−16281. (21) Wheeler, D. R.; Rowley, R. R. Mol. Phys. 1998, 94, 555−564. (22) Yeh, Y. L.; Zhang, C.; Held, H.; Mebel, A. M.; Wei, X.; Lin, S. H.; Shen, Y. R. J. Chem. Phys. 2001, 114, 1837−1843. (23) Pártay, L.; Jedlovszky, P.; Horvai, G. J. Phys. Chem. B 2005, 109, 12014−12019. (24) Picaud, S.; Hoang, P. N. M. J. Chem. Phys. 2000, 112, 9898− 9908. (25) Hantal, G.; Jedlovszky, P.; Hoang, P. N. M.; Picaud, S. Phys. Chem. Chem. Phys. 2008, 10, 6369−6380. (26) Schuler, L. D.; Daura, X.; van Gunsteren, W. F. J. Comput. Chem. 2001, 22, 1205−1218. (27) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. J. Phys. Chem. 1987, 91, 6269−6271. (28) Mezei, M.; Beveridge, D. L. Ann. Acad. Sci. N.Y. 1986, 482, 1− 23. (29) Leach, A. R. Molecular Modelling; Longman: Singapore, 1996. (30) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935.

Figure 4. Dependence of the D parameter (see eq 19) on the acetone mole fraction, as obtained in mixtures of TIP5P-E water and PAC acetone by direct double derivation of the free energy of mixing data (filled circles), as well as by fitting the Margules equation (eq 20) to these data and performing the double derivation analytically (solid curve).

different model combinations revealed that lower Amix values are obtained with an acetone model of higher, or with a water model of lower, dipole moment. The free energy cost of mixing of any model pair in any ratio turned out to be very small, being below 0.5RT, i.e., the mean kinetic energy of the molecules along one degree of freedom. All these findings stress the difficulty of accurately describing the thermodynamic changes accompanying the mixing of acetone and water and point out the key role of the description of the electrostatic interactions in this respect. Unlike the other acetone models, the recently proposed model of Pereyra, Asar, and Carignano10 turned out to be indeed miscible with at least the TIP5P-E model36 of water in any ratio. Moreover, this model combination reproduced the change of several thermodynamic properties (i.e., energy, free energy, and TS) rather accurately, within 0.8 kJ/mol at any mole ratio. Furthermore, the agreement between the simulation results and experimental data, concerning at least the energy of mixing, can be even further improved by treating the acetone molecules as flexible.10 Moreover, comparison of the free energy of mixing of TIP5P-E water with other water models, considering either the AUA4 or the TraPPE model of acetone (see Figure 2), suggests that the PAC acetone model is likely to be fully miscible also with other water models, at least with SPC and TIP4P. The PAC model, parametrized to reproduce the enthalpy of mixing with TIP5P-E water,10 accounts for the polarization of the acetone molecules due to their strongly polar water neighbors in an average way, as the fractional charges of this model scale with the mole ratio of acetone−water mixtures. Although this treatment of the acetone polarization evidently restricts the use of this model only in acetone−water binary mixtures, its success in reproducing the mixing behavior of acetone and water in both a qualitative and quantitative sense clearly stresses the importance of taking the effect of molecular polarization into account in the modeling. The present results also point out the urgent need for the development of a new, general purpose acetone model that is also able to reproduce full miscibility with water and strongly suggest that such a model needs to be polarizable. 5983

dx.doi.org/10.1021/jp302629r | J. Phys. Chem. B 2012, 116, 5977−5984

The Journal of Physical Chemistry B

Article

(31) Kang, M.; Perera, A.; Smith, P. E. J. Chem. Phys. 2009, 131, 157101−1−2. (32) Jedlovszky, P.; Idrissi, A.; Jancsó, G. J. Chem. Phys. 2009, 131, 157102−1−2. (33) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed; Reidel: Dordrecht, 1981; pp 331−342. (34) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505−1− 12. (35) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910−8922. (36) Rick, S. W. J. Chem. Phys. 2004, 120, 6085−6093. (37) Foloppe, N.; MacKerrel, A. D., Jr. J. Comput. Chem. 2000, 21, 86−104. (38) Martin, M. G.; Biddy, M. J. Fluid Phase Equilib. 2005, 236, 53− 57. (39) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (40) Darvas, M.; Jedlovszky, P.; Jancsó, G. J. Phys. Chem. B 2009, 113, 7615−7620. (41) Boje, L.; Hvidt, A. J. Chem. Thermodyn. 1971, 3, 663−673. (42) Mezei, M. MMC Program at URL: http://inka.mssm.edu/ ∼mezei/mmc. (43) Barker, J. A.; Watts, R. O. Mol. Phys. 1973, 26, 789−792. (44) Neumann, M. J. Chem. Phys. 1985, 82, 5663−5672. (45) Smith, J. M.; Van Ness, H. C.; Abbott, M. M. Introduction to Chemical Engineering Thermodynamics, 6th ed.; McGraw Hill: Boston, 2001.

5984

dx.doi.org/10.1021/jp302629r | J. Phys. Chem. B 2012, 116, 5977−5984