- Email: [email protected]
A comparative Molecular Dynamics study of water–methanol and acetone–methanol mixtures
Journal of Molecular Liquids 159 (2011) 52–59
Contents lists available at ScienceDirect
Journal of Molecular Liquids j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / m o l l i q
A comparative Molecular Dynamics study of water–methanol and acetone–methanol mixtures Aurélien Perera a,⁎, Larisa Zoranić a,b, Franjo Sokolić b, Redha Mazighi a a b
Laboratoire de Physique Théorique de la Matière Condensée (UMR CNRS 7600), Université Pierre et Marie Curie, 4 Place Jussieu, F75252, Paris cedex 05, France Department of Physics, Faculty of Sciences, University of Split, Nikole Tesle 12, 21000, Split, Croatia
a r t i c l e
i n f o
Article history: Received 3 February 2010 Received in revised form 20 April 2010 Accepted 7 May 2010 Available online 13 May 2010 Keywords: Aqueous mixtures Molecular liquids Simulations Statistical mechanics
a b s t r a c t Two binary mixtures with methanol as common component, namely water–methanol and methanol– acetone are studied by Molecular Dynamics simulations. Thermodynamical properties such as enthalpies, excess enthalpies, volumes and excess volumes are compared. Structural properties are studied through the various site–site pair correlation functions and the associated Kirkwood–Buff integrals. While the thermodynamical properties are relatively well calculated, we show that structural properties are affected by two severe problems. The first is an inherent property of the asymptote of the correlation function in finite systems, which affects their integrals, the second is the micro-heterogeneous structure of these hydrogen bonding mixtures, which affects the medium-to-long-range part of the distribution functions. These two problems conspire to make computer simulation unreliable to study this precise part of the structural properties of such mixtures. We propose pathways to correct this situation and demonstrate how it serves to considerably improve the calculated Kirkwood–Buff Integrals. Finally, we compare the two mixtures relative to their tendencies to form local heterogeneity. The analysis demonstrates that the microstructure of neat methanol is better preserved in acetone than in water. © 2010 Elsevier B.V. All rights reserved.
1. Introduction It is now relatively well documented and accepted that mixtures of associating liquids, such as aqueous mixtures, are difficult to simulate, and in particular, need excessively large simulation times and system sizes [1–5], which is a real paradox, mainly in view of the small size of these molecules. The reason for this unusual behaviour is that molecules preferentially self-associate through the strong hydrogen bonding interaction, and therefore, in addition to the molecular dynamical scale, they introduce a second dynamic scale — that of the clusters that are formed. As a result, one needs large system sizes, mainly to probe the statistics of the various clusters that are formed, and because of the slower motion of the clusters dynamics, one also needs excessively long run times. Paradoxically, this problem is less serious for the neat systems, such as water or alcohols, where, nevertheless, considerable micro-structure have been shown to exist [6–9]. The reason for this difference may be tracked back to the fact that clusters in neat systems have either a well defined geometry, such as those in neat alcohols [8], or are too elusive to be detected by ordinary statistics, such as water or amides [10]. In contrast, mixtures of such liquids tend to have some sort of self-partitioning [12,13], quite visible in snapshots, which confers these systems a rigidity that ⁎ Corresponding author. E-mail address: [email protected] (A. Perera). 0167-7322/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2010.05.006
defies simulation techniques. Nanostructuring is not unique to aqueous mixtures and can be observed in computer simulations of various ionic liquids as well [11], where hydrophobic chains selfsegregate from the ionic headgroups. There is an extensive computer simulation studies of aqueous methanol mixtures by several authors [3,12–17], and this system is now recognized to represent a typical example of well mixed system that exhibits strong micro-heterogeneity (MH). We revisit here this mixture, that we will consider as a reference system, in order to compare with the acetone–methanol mixture, which is supposed to be less structured, from the simple argument that methanol has only one donor and one acceptor bonding sites, as opposed to water which has two of each. In an earlier work [18], we have studied the acetone– water mixture and shown that most existing models lead to demixing schemes. In view of this unsatisfactory status, we study here the acetone–methanol mixture, which is well behaved with respect to the demixing issue. It is noteworthy that even this simple mixture is not satisfactorily simulated with usual force field generation techniques and specific issues have been discussed in recent works [15–17]. There are less studies of the acetone–methanol mixtures, to our knowledge. Computer simulation studies [19,20] of systems of less than 600 molecules report the existence of local heterogeneity and methanol chain clustering. This appears to be confirmed by various spectroscopic techniques [21–23], although the extent of segregation differs according to different techniques [22,23].
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
The main focus of this paper is how micro-heterogeneity influences the calculation of the structural functions and how to deal with this problem, and when this is possible to deal with it efficiently. For this reason, it is important to examine mixtures with expected different behaviour and see how they are affected by their respective micro-heterogeneity. 2. Simulation and theoretical details 2.1. Simulations In this work, we have considered OPLS force field models for both acetone [24] and methanol [25], and SPC/E model for water [26]. We emphasize again that neither SPC/E nor TIP4P water could mix satisfactorily with OPLS acetone [18]. The problem of detecting a true demixing within finite size computer simulation is not easily resolved when strong micro-heterogeneity is present. This important issue has not been satisfactorily addressed at present. The problematic is the following. Phase separation is induced by strong concentration fluctuations. These can be tackled properly within computer simulations by scaling techniques [27]. However, when local heterogeneity is present, it becomes difficult to disentangle this feature from genuine concentration fluctuations. For both mixtures, we have performed Molecular Dynamics simulation in the isobaric ensemble, for the ambient conditions of temperature and pressure. Berendsen thermostat and barostat where applied, with relaxation times of 0.1 ps and 0.5 ps, respectively. The time step for the integration of motion was fixed at 2.10− 15 s in all the calculations, which was found satisfactory in order to reproduce the rotational motion of water. In all our simulations total of N = 2048 molecules is considered for each mixture. This size allows us to reach intermolecular distance between 2.5 nm and 3 nm, which is about 8 to 10 times the molecular size, and which should be, in principle, sufficient to reach proper asymptotic bulk behaviour for the pair correlation functions. This is not often true, as we discuss in the results section. For each mixtures, the mole fraction of one of the component was varied by steps of 0.1, from zero to one. Each such system was equilibrated for 100–120 ps. This time window is sufficient to establish equilibrium values for the energies and volumes. Statistics were performed in batches of run length varying in the range 64– 256 ps, for about 10–20 batches. The accumulated statistics cover a range from 1 to 3 ns. This is still about twice smaller than what is claimed to be optimal in recent works [1]. The DLPOLY2 code [28] was used to perform these simulations. As pointed out in our earlier work [5], this code does not reproduce the expected asymptotic behaviour of the site–site correlation functions, and a multiplicative correction factor N/(N − 1) needs to be applied. This problem is independent of the issue of the asymptotical behaviour of these functions inherent to finite size, that we now expose. 2.2. Handling the asymptotic behaviour of pair correlation functions in mixtures Although computer simulations are performed with periodic boundary conditions, they are nevertheless not handling infinite size system. It is a well known result of statistical mechanics [29] that the asymptotical behaviour of the radial distribution function of an finite system, for example the canonical or isobaric ensemble, goes to a limit which is 1 − 1/N rather than exactly one as it should be for an open system (grand canonical ensemble). In the thermodynamic limit, which is when N (the number of particles) and V (the volume) go to infinity, then the equivalence of all these ensembles is automatically achieved [29]. But, in practice, the problem of the asymptote not being unity is a cumbersome one if one wants to compute the Kirkwood–Buff integrals (KBI), which are integrals over the pair distribution functions. Surprisingly, this important point has not been pointed out in many
53
earlier works concerning the evaluation of the KBI. But there is even more. In many textbooks, the limit that is proposed, which is 1 − 1/N, is in fact that of the ideal gas system. In presence of interactions and correlations, the proper result has been known for decades from the earlier works [30,31] of Lebowitz and Percus in particular, and it is 1 − 1/(NκT⁎), where κT⁎ = ( ∂ βP / ∂ ρ)T is the reduced isothermal compressibility, where β = 1/kBT is the Boltzmann factor, P the pressure and ρ the density. For the ideal gas, κT⁎ = 1, and the often quoted limit is recovered. Otherwise, there is a clear asymptote to below unity, that can create problems if the KBI are evaluated from this RDF. The extension of this expression for mixtures [32], in terms of the pair correlation between the pair of species i, j, reads: lim g ð1; 2Þ r→∞ ij
1 1 ∂ρ = 1− qffiffiffiffiffiffiffiffiffiffi ð i ÞTNk Ni Nj ρ ∂βμ j
ð1Þ
where μi is the chemical potential of species i. The evaluation of this derivative in practical calculations is a tedious task, since a chemical potential is needed, which is already a very noisy quantity to obtain from simulations, hence its derivative would be entitled to even larger errors. It is important to realize that, for the case of mixtures, the 1/N factor does not necessarily imply that the asymptotical value will be small if N is large. Indeed, for small concentration of one of the species, the corresponding factor NiNj can be arbitrarily small. We illustrate this in the results section. The conclusion of this fact, is that the asymptote of the pair distribution function can only be revealed by looking closely at the calculated functions. The expression above holds equally, as such, for the site–site correlation functions, since all site– site functions integrate out to give the same KBI as that of the centerof-mass correlation function [29]. This result is a consequence of the invariance of this quantity on the choice of the molecular center [29]. We propose a simple way to cure this problem artificially. We first spot the asymptotical behaviour of the function by direct visualization. Then, we add the following correcting function to the correlation function: " correct gij ðrÞ
=
raw gij ðrÞ
δij 1+ 2
r−κij 1 + tanh αij
!# ð2Þ
where the coefficients δij, κij and αij depend on the considered site–site correlation function. Usually, we take κij to be twice the distance at which gij(r) starts to be different from zero (an effective diameter). δij corresponds to the shift observed in the incorrect asymptote value aij with δij = 1 − aij, and αij controls the smoothness of the passage from aij to 1. We take this parameter to be a constant factor about twice the molecular diameter. We have verified that this value could be varied from the molecular diameter to about 3 times this distance without noticeable changes in the integrals. The reason of this invariance is due to the fact that the amplitude of the change is of the order of 0.1%. However small this value could be, it nevertheless leads to noticeable deviation of the KBI as will be illustrated in the next section. The form of the function above avoids perturbing the original function near the first few neighbours, and in particular in the core region, which is an absolute requirement. Finally, we note that, for neat systems, this tail correction can be so small, usually when system sizes larger than N = 1000 are considered, that the calculated RDF tend to unity within the numerical error [33]. This is particularly the case when the compressibility of the neat liquid is very small, as it is often the case with most liquids. 3. Results We present thermodynamic properties separately from the structural ones principally because the first ones do not seem to be affected the same way as structural ones. Properties such as enthalpies
54
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
and volumes converge much faster towards their equilibrium values and need only modest system sizes, such as N = 300–500 molecules. This is in sharp contrast with the proper computation of the radial distribution functions, which requires larger systems and long run times. This important fact indicates that the strong local inhomogeneity has a definitive impact on the evaluation of different properties. This is in contrast with the evaluation of the same properties in simple Lennard–Jones liquids and mixtures, where the statistics of both quantities require the same spatial and temporal scales, except near phase separation boundaries. 3.1. Thermodynamic properties We present the enthalpies and the volumes, as well as their excess quantities defined for the quantity A(x) as Aex(x)=A(x) −A1 −x(A2 − A1), as a function of x, the mole fraction of species labeled 2, and A1 =A(0) and A2 =A(1) are the neat liquid values of species 1 and 2, as obtained for x=0 and x=1, respectively. We compare our results with experiments and other works for the case of aqueous methanol solution. The reported enthalpies are in fact configurational energies, since we maintain atmospheric pressure, the PV term is negligible [34]. Fig. 1 shows the densities (top panel), volumes (middle panel) and excess volumes (lower panel) for the aqueous methanol mixture (labeled from now WM) while Fig. 2 shows the same quantities for the acetone–methanol mixtures (labeled AM). There is an obvious difference in the excess quantities: that of the AM mixture is nearly zero, indicating very little volume change when mixing the two components. Since mixing two hard spheres liquids leads to a noticeable negative excess volume [35], this means that the flat excess volumes of AM mixture indicate an increasingly repulsive excess volume as acetone is added into the mixture. Since hard spheres mix homogeneously, none of the tendencies observed here demonstrate any a priori peculiar behaviour related to the formation of micro-heterogeneity. We recall that
Fig. 1. The densities (top panel), volumes (middle panel) and excess volumes (lower panel) for aqueous methanol , as functions of the methanol mole fraction. The red curve is experimental data from Ref. [36] and excess volumes from Ref. [37]. Black symbols are for the present MD data, green symbols from simulations of Ref [3] and magenta symbols from Ref. [14]. In the lower panel blue symbols are from Ref. [17].
Fig. 2. The densities (top panel), volumes (middle panel) and excess volumes (lower panel) for acetone–methanol mixture, as functions of the acetone mole fraction.
the aqueous acetone excess volumes were found negative for many model combinations while most showed tendency for demixing [18]. Fig. 3 shows the energies and excess energies for the WM mixture and Fig. 4 for the AM mixture. The principal difference between the two mixtures is that the excess for WM are negative while they are positive for the AM mixture. This may be thought to be consistent with the behaviour found for the excess volumes: the reason for near flat excess
Fig. 3. The energies (top panel) and excess energies (lower panel) for aqueous methanol, as functions of the methanol mole fraction. Line and symbols as in Fig. 1. The red line is experimental data from Ref. [38].
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
55
MH problem, and then by computing the KBI, which are thermodynamic quantities that relate directly to probing through measurable quantities such as partial molar volumes and chemical potentials. 3.2.1. Snapshots Fig. 5 shows a snapshot of the aqueous methanol mixture for 30% methanol content, where the water(left) and methanol(right) molecules are shown separately for the same angle of shot. It is visually obvious that water molecules have a Swiss-cheese like distribution, noticeable from the fact that distant molecules are seen as smaller than those closer to the viewer. This is particularly obvious when the shots are made to rotate when visualizing them. In contrast, the methanol molecules appear more homogeneously distributed. This means that some of the water clusters contain embedded methanol molecules. Similar snapshot-based conclusions about micro-segregation were reported in other works [12] and for other mixtures [1,10,18]. Fig. 6 shows snapshots of the acetone–methanol mixture for 30% acetone content. The micro-heterogeneous distribution of the methanol molecules (left) is quite noticeable, in contrast to that of the acetone molecules (right) which appears more homogeneous. These observations indicate that the common species, methanol, plays a different role in different environment. In particular, these shots seem to imply that it is always the more strongly associating molecule that drives the micro-segregation. Fig. 4. The energies (top panel) and excess energies (lower panel) for the acetone– methanol mixture, as functions of the acetone mole fraction.
volumes of the AM mixture can come from the overall repulsiveness witnessed through the excess energies. We recall that the excess energies for the aqueous acetone mixtures were found positive for nearly all mixtures [18]. In contrast, however, the AM mixture appears well mixed, despite the inherent micro-heterogeneity. These various trends do not allow to guess the micro-heterogeneous status of these mixtures, similarly to that found for the volumes. We note that, for the aqueous methanol mixture, most calculated energies are in relatively good agreement with the experimental ones, except for the fact that all calculated excess energies have a minimum at x = 0.5 while experimental ones are positioned at x ≈ 0.3. This strange fact has not been yet satisfactorily interpreted. 3.2. Structural properties The study of these properties can be made in three steps. First by showing evidence of MH through the snapshots, next by studying the correlation functions themselves, which is a statistical approach to the
3.2.2. Distribution functions Fig. 7 shows selected site–site radial distribution functions (ssRDF) to illustrate the structure of the aqueous methanol mixture for methanol mole fractions of 0.20, 0.50 and 0.80. These plots indicate that the sensitivity of the water–water and meth–meth correlations, to the variation of methanol mole fraction, is quite different. Water– water correlations (in blue) show a dramatic increase of the first few neighbours structure, which indicates that water forms tight clusters in methanol rich environment. In comparison, the various methanol– methanol correlations seem relatively insensitive to changes in water concentration. In particular, the methanol oxygen–oxygen (O–O) correlations show the inverse trend to that of water: they increase with methanol concentration while the water–water correlations increase with decrease of water concentration. Since both species are hydrogen bonding species, one would expect a symmetric behaviour, and the observed asymmetry is rather striking. The interpretation of this fact corroborates the snapshots of Fig. 5. Water tends to microsegregate more than methanol, resulting in an increased correlation when it is scarce in methanol, as it tends to form small tight clusters. Conversely, methanol is more homogeneously distributed in water, as
Fig. 5. Snapshots of the aqueous methanol mixture at 30% methanol content. Water molecules are shown in the left and methanol molecules at the right. In both shots the missing molecules are shown in skeleton. The color code is white for H atoms, red for O atoms and blue for the methyl group.
56
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
Fig. 6. Snapshot of the acetone–methanol mixture at 30% acetone content. Methanol molecules on the left and acetone molecules on the right. The color codes are the same as in Fig. 5.
can be seen from the little variation of the methyl group correlations. There must be some micro-segregation of methanol at low content (20%), since its distribution is nearly that of 80% methanol content, but one cannot really speak of methanol cluster the same way water is clustered in 80% methanol, because of the marked differences between the first peaks in the respective cases. By comparing the structure between the different meth–meth correlations we can guess the way methanol molecules orient in various conditions. In neat methanol, we have previously compiled the structure of different methanol clusters [9,32] and noted that the methyl (M) sites were homogeneously distributed while the O sites were connected in
Fig. 7. Selected site–site distribution functions for the aqueous methanol mixture, at methanol mole fraction x = 0.20, 0.50 and 0.80. The oxygen–oxygen (OO) are shown in blue for water, black for methanol and magenta for water–methanol. The red curve is for methanol MM correlations and the green curve for water–Methanol OO cross correlations.
various chain-like patterns. Assuming this is still the case in 80% meth content, we see that the decrease of the O–O correlations with methanol content indicates that methanol O sites would preferentially bond with water. The increasing cross correlations between the oxygens of the two species (in magenta), with methanol decrease, confirm this hypothesis. In contrast, the relative insensitivity of the M–M correlations (in red) indicate that the methyl sites are homogeneously distributed within the methanol rich areas in all circumstances. These various trends confirm the local micro-heterogeneous state of the mixture. Fig. 8 shows selected ssRDF for the acetone–methanol mixture. Here, only methanol is self-bonding. Therefore, it is not surprising to see methanol O–O correlations (in magenta) behave similarly to water–water correlations in Fig. 7. Indeed, the first neighbour correlation increases dramatically with acetone concentration, due to strong clusterisation of methanol in acetone rich environment. In comparison, we see that the acetone O–O correlations (in blue) are insensitive to changes in concentration, indicating that acetone is correlated to itself much the same way in concentrated and diluted situations. This is in favour of micro-segregation of the two species. The methanol M–M correlations (in red) increase when methanol decreases, although less dramatically than the O–O correlations. This indicates that the M sites are not so randomly distributed (as they appeared to be in aqueous methanol), enforcing the picture that the methanol clusters should be quite similar to that of the neat methanol (chain-like). The O–O cross correlations (in green) indicate that little Hbonding occurs between the two species. These two figures illustrate the difference in micro-heterogeneity in the two types of mixtures sharing a common Hbonding molecule. From the detailed analysis above, the major difference is in the cluster structure of methanol in the two different environments. The microstructure of neat methanol clusters appears to be more preserved in acetone than in water. We turn now to the asymptotical behaviour of these distribution functions. Fig. 9 serves to illustrate the fact that the asymptote of the correlation functions do not go to unity, but behave rather like in Eq. (1). For this purpose we pick the acetone–methanol mixture at 20% acetone content, and plot the following two methanol–methanol ssRDF, oxygen–oxygen and methyl–methyl. The original functions are shown in dashed lines and those corrected for the asymptote in full lines. The top panel shows these 2 functions and the differences with the corrected ones cannot be distinguished on this scale. The middle panel shows a zoom on the tail part of the ssRDF. The curve in red shows the correcting function suggested in Eq. (2), and the green curve shows the current asymptote, which we estimate at a = 0.9983 from the plot, which is obviously different from unity. It is important to note that both ssRDF tend to the same asymptote. The curves
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
Fig. 8. Selected site–site distribution functions for the acetone–methanol mixture, at acetone mole fraction x = 0.20, 0.50 and 0.80. The magenta curve is for meth–meth O–O correlations, the red curve for meth–meth M–M correlations, the blue curve for acetone–acetone O–O correlations and the green curve for acetone–methanol O–O cross correlations.
corrected for the asymptote are shown in full lines, with the following parameters in Eq. (2), κ = 5Å and α = 1 Å. The lower panel shows the corresponding running KBIs (RKBI) defined as [10] r
2
Gij ðrÞ = 4π∫0 dt t ðgij ðrÞ−1Þ
ð3Þ
where gij(r) is the ssRDF between sites i and j. The asymptotic behaviour of this latter function, when r → ∞ is precisely the KBI. There is a major difference between corrected and uncorrected KBIs, so the correction is really needed. The downwards curvature of the uncorrected curve (in dashes) comes simply from the integrated r3 behaviour in Eq. (4), that arises if the asymptote is not strictly unity. It is important to note that this correction is not due to the existence of MH, but to the finite size effect, and can also be seen in simple liquids [31] and their mixtures [32]. However, the existence of MH can mask this behaviour, depending on the extension of correlations which are affected by the fact that this micro-heterogeneity extends up to the half box size, which is the usual cutoff-distance for the computation of the correlations. Finally, we show an example of the behaviour of the correlation function where the existence of MH does not allow a proper settling of the asymptotical behaviour of the RDF and hence of the corresponding RBKIs. This is illustrated in Fig. 10, where we show selected ssRDF for the equimolar aqueous methanol mixture. Each curve presents the outcome of a batch of statistics for run length of 64 ps. The average that we consider acceptable is shown in thicker black line. It is clear that there is no convergence of the behaviour to any proper definite value. This type of problem is recurrent in the simulation of any aqueous mixture, and is never reported in the literature. It should also occur for other microsegregating functions, such as for aqueous ionic-liquid mixture for
57
Fig. 9. (Top panel) The oxygen–oxygen (OO in blue) and methyl–methyl group (MM in magenta) site–site correlation functions for the acetone–methanol mixture with 20% acetone. (Middle panel) Zoom on the asymptotical behaviour of the OO and MM correlation of the top panel. The correcting function (Eq. (2)) is shown in red (parameters in the text). The current asymptote is shown in green. The corrected functions are shown in full lines and the uncorrected ones in dashes. (Lower panel) The RKBIs corresponding to the functions above, with the same color conventions.
example. These small differences that are amplified in the evaluation of the RKBI are unnoticeable in the first few peaks of the full RDF (upper panel). This could be the reason why these problems are never noticed in the published literature, when showing the RDF under the 12 Å range. We wish to emphasize that correlations at the first and second peak are relatively well stabilized after few 100 ps, and it is only the behaviour beyond that is subject to the anomalous behaviour reported here. Hence, results reported for smaller system sizes should give a reliable account of the correlations of the first two neighbour shells. But one should keep in mind that this fact masks the inherent pathological behaviour for higher neighbour shells, that is ultimately due to a real physical effect, that of the micro-heterogeneity. 3.2.3. Kirkwood–Buff integrals In Fig. 11, we report the KBI of the two mixtures, together with the experimental results for the aqueous methanol mixture. We are not aware of any experimental results for the acet–meth mixture. In place, we have plotted an estimate of the experimental KBIs, evaluated from the following expressions [32]: G12 ðxÞ = −
V1 V2 VðxÞDðxÞ
Gii ðxÞ = G12 ðxÞ +
Vj 1 −VðxÞ xi DðxÞ
ð4Þ
where x = x2 is the mole fraction of component 2 (methanol in water– meth mixture and acetone in acet–meth mixture), V1 and V2 are the pure component volumes, V(x) = V1 + x(V2 − V1) is the approximate
58
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
Fig. 10. (Top panel) The water–water radial distribution function of aqueous methanol equimolar mixture for different statistics (undistinguishable on this scale). (Middle panel) Zoom on the asymptotical behaviour of the RDFs above, each curve represents a statistics over 64 ps. (Lower panel) The RKBIs corresponding to the functions above, with the total average shown in the thicker black curve.
volume of the mixture and D(x) = 1 − x(1 − x)d is related to the chemical potential derivative in a crude way, and i = 1,2 with j = 2,1 in the second formula. These expressions, which neglect the compressibility of the mixture as well as specificity of the partial molar volumes and chemical potentials, reproduce quite well the experimental KBIs of the aqueous methanol mixture. It is often found that modifying the definition of D(x) can allow a very good fit of the KBIs [39]. The volume for the neat water is taken as V = 18 cm3/mol, for methanol V = 40 cm3/mol and for acetone V = 78 cm3/mol. We use d = 1.20 for aqueous methanol and d = 2.10 for the acetone–methanol mixture. The MD results are obtained by correcting the asymptotic result whenever this is possible. As noted above, this is possible when the asymptote is not clouded by bad statistics due to MH. We find that this was invariably the case when considering self-correlations of the majority species, as illustrated in Fig. 9. Conversely, self-correlation of the minority species was always difficult to correct since the statistics of the clustering strongly affect the stabilisation of the asymptote. In such bad cases, the asymptotical behaviour is chosen by averaging the different values (see Fig. 10). The reported results indicate that the KBI evaluated after correction are nevertheless in very good agreement with the experimental results for the aqueous methanol mixture and with the fitting expressions of Eq. (4) above. The first thing that we notice is that the KBIs are relatively well defined for the solute (component 2) in almost the full solute mole fractions, while for the solvent this is true only in the solvent rich region. In general, the component poor region leads to poor estimate of the corresponding KBIs. This is clearly due to the large clustering of such species in the adverse rich component, which in turn affects the behaviour of the tail of the ssRDF. We note that methanol in the two mixtures has a similar behaviour when it is in low concentration: the KBI are noisy and out of bounds. This is again due to clustering and therefore bad definition of the asymptote. We note again that methanol, that is the common species, behaves very differently when it is in water than when it is in acetone. The striking similarity of the behaviour of methanol in aqueous methanol with that of acetone in acetone–methanol, as well as that of methanol in acetone–methanol with that of water in aqueous methanol, emphasizes the asymmetrical role of this substance in different environments. These particularities are due to the Hbonding nature of the species. The next obvious remark is that there is no direct signature of the differences in the micro-heterogeneous behaviour of these two mixtures that can be detected from simple visualization of the KBIs. Therefore, the ssRDF remain the best observables to appreciate the MH of a mixture, albeit the problems encountered in evaluating them. 4. Discussion
Fig. 11. The KBI for the aqueous–methanol mixture versus methanol mole fraction (top panel) , and acetone–methanol mixtures versus acetone mole fraction(lower panel). MD results are shown in filled dots (top panel: blue for water–water, red for Meth– Meth and black for Water–Meth; bottom panel: blue for Meth–Meth, red for Acet–Acet and black for Meth–Acet). Experimental results from Ref. [39] are shown in the top panel in open squares with color convention similar to the MD symbols. The fit from Eq. (4) are shown in lines with color convention similar to the MD symbols.
The fact that the simple asymptote correction gives appreciably good KBI indicates the importance of such correction. It is also obvious that rather large system sizes must be studied, as can be clearly seen from the zoom of Figs. 9 and 10: one needs to reach the part where oscillations due to packing structure die out, and the bulk regime starts to show up. Generally, this happens roughly after 5 molecular diameters, although no general rules can be easily specified. The paradox we pointed out earlier, namely that one needs such large systems and therefore long statistics and runs for such small molecules mixtures, is truly disconcerting one. We would like to stress that these calculations were conducted with unmodified versions of the OPLS force fields of the methanol and acetone molecules. The fact that a good expectation of the KBIs is obtained in this work, enforces the idea that it may not be really necessary to modify existing force fields, such as the OPLS force field for example, in order to achieve this goal. Such alternate routes were indeed taken recently with success by some authors [1–4]. The present study
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
shows that the painstaking task of redesigning force fields could be avoided by explicitly taking into account inherent corrections due to finite size effects. Following this route, it seems that such type of corrections could be extended to account for the spurious contributions due to the local heterogeneity. We would like to stress, however, that in cases where spurious demixing is enforced by improper force fields [18], it would be rather difficult to use such tricks to improve results on KBI, for obvious physical reasons. Ultimately, the pathways followed in this work amount to distinguish between contributions due to concentration fluctuations and those due to local micro-heterogeneity [5]. As stated in the Introduction, the first contributions control phase separation while the second are responsible for the difficulties encountered in simulations of these particular systems. It is at present quite difficult to predict how these two contributions stem from the force field parametrisation. Recently, some authors [16,17] have raised the subtle issue of reproducing the experimentally observed minimum in the excess partial molar volume of the aqueous alcohol mixtures. It seems that none of existing force fields are able to reproduce such minimum. Another, perhaps related issue, is the shift in minimum in the experimental excess enthalpy of aqueous methanol around x = 0.28, instead of 0.5 for all force field models (see lower panel of Fig. 3). Ultimately, a convincing computer simulation study would have to be able to account for the subtle local structural orders that can give rise to these still unexplained features. 5. Conclusion This paper focuses on the micro-heterogeneous structure of different types of hydrogen bonding mixtures and the ability to report on this property through MD simulations. For this purpose, two mixtures have been studied with a common component, methanol. The main conclusion of this work is that these micro-segregating mixtures can be studied by MD simulation in a reliable way if one focuses on properties such as enthalpies and volumes, as well as correlation functions when these are limited to the first and second neighbours shells. Generally speaking, these correlations tend to die rather quickly after the second and third neighbours range. This is why most reported results in the literature are never concerned with this problematic. However, this apparent success should not mask the real underlying problem caused by the fact that these mixtures are strongly inhomogeneous, and that such inhomogeneity affects undeniably long range properties of these mixtures in an uncontrollable way. This is due to two facts: first the inherent problem of the asymptote, for which we provide here a remedy, and second the extent of the local clustering that can affect the convergence of the value of the asymptote. Increasing the size of the system may appear as a remedy, but it comes with additional
59
problems such as the need for longer runs, since one should now also take into account the slow dynamics of the clustered parts of the simulated sample. The present calculations show that methanol behaves differently in the aqueous and acetone environments, mainly with respect to clustering. It clusters less in water than in acetone. This fact alone indicates that micro-heterogeneity is not an intuitively predictable property of these mixtures.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]
R. Chitra, P.E. Smith, J. Phys. Chem. B104 (2000) 5854. R. Chitra, P.E. Smith, J. Chem. Phys. 115 (2001) 5521. S. Weerasinge, P.E. Smith, J. Chem. Phys. 118 (2003) 10663. M.E. Lee, N.F. van der Vegt, J. Chem. Phys. 122 (2005) 114509. L. Zoranic, R. Mazighi, F. Sokolic, A. Perera, J. Chem. Phys. 130 (2009) 124315. J.-H. Guo, Y. Luo, A. Augustsson, S. Kashtanov, J.-E. Rubensson, D.K. Shuh, H. Agren, J. Nordgren, Phys. Rev. Lett. 91 (2003) 157401. R. Ludwig, ChemPhysChem 6 (2005) 1369. L. Zoranic, F. Sokolic, A. Perera, J. Chem. Phys. 127 (2007) 024502. A. Perera, F. Sokolic, L. Zoranic, Phys. Rev. E75 (2007) 060502-(R). L. Zoranic, R. Mazighi, F. Sokolic, A. Perera, J. Phys. Chem. C111 (2007) 15586. J.N.A. Canongia Lopes, A.A.H. Padua, J. Phys. Chem. B110 (2006) 3330. S. Dixit, J. Crain, W.C. Poon, J.L. Finney, A.K. Soper, Nature 416 (2002) 829. L. Dougan, S.P. Bates, R. Hargreaves, J.P. Fox, J. Crain, J.L. Finney, V. Reat, A.K. Soper, J. Chem. Phys. 121 (2004) 6456. A.M. Ferrario, M.M. Haughney, I.R. McDonal, M.L. Klein, J. Chem. Phys. 93 (1990) 5156. E.J.W. Wensik, A.C. Hoffmann, P.J. van Maaren, D. van der Spoel, J. Chem. Phys. 119 (2005) 7308 52003. L. Vlček, I. Nezbeda, J. Mol. Liq. 131–132 (2007) 158. Diego González-Salgado, I. Nezbeda, Fluid Phase Equilib. 240 (2006) 161. A. Perera, F. Sokolic, J. Chem. Phys. 121 (22) (2004) 11272. D.S. Venables, C.A. Schmuttenmaer, J. Chem. Phys. 113 (2000) 3249. G. Kamath, G. Georgiev, J.J. Potoff, J. Phys. Chem. B109 (2005) 19463. D.S. Venables, C.A. Schmuttenmaer, J. Chem. Phys. 113 (2000) 3243. M. Musso, M.G. Giorgini, H. Torii, J. Mol. Liq. 147 (2009) 37. J.-J. Max, C. Chapados, J. Chem. Phys. 122 (2005). W.L. Jorgensen, J.M. Briggs, M.L. Contreras, J. Phys. Chem. 94 (1990) 1683. W.L. Jorgensen, J. Phys. Chem. 90 (1986) 1276. J.C. Berendsen, J.P.M. Postma, W.F. Von Gusteren, J. Hermans, in: B. Pullman (Ed.), Intermolecular Forces, Reidel, Dortrecht, 1981. N.G. Wilding, Phys. Rev. E. 52 (1995) 602. DL_POLY package from CCLRC Daresbury Laboratory www.ccp5.ac.uk/DL_POLY/. J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986. J.L. Lebowitz, J.K. Percus, Phys. Rev. 122 (1961) 1675. J. Kolafa, S. Labik, A. Malijevsky, Mol. Phys. 100 (2002) 2629. A. Perera, L. Zoranic, F. Sokolic and R. Mazighi (in preparation). L. Zoranic, F. Sokolic, A. Perera, J. Mol. Liq. 136 (2007) 199. F. Sokolic, A. Idrissi, A. Perera, J. Chem. Phys. 116 (2002) 1636. J.S. Rowlinson, F.L. Swinton, Liquids and Liquid Mixtures, Butterworths, 1982. Landolt-Börnstein, Neue Serie Group IV , vol.1 part B, SPringer Verlag, 1977. Y. Koga, Can. J. Chem. 77 (1999) 2039. L. Benjamin, G.C. Benson, J. Am. Chem. Soc. 67 (1963) 858. A. Perera, F. Sokolic, L. Almasy, Y. Koga, J. Chem. Phys. 124 (2006) 124515.
Contents lists available at ScienceDirect
Journal of Molecular Liquids j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / m o l l i q
A comparative Molecular Dynamics study of water–methanol and acetone–methanol mixtures Aurélien Perera a,⁎, Larisa Zoranić a,b, Franjo Sokolić b, Redha Mazighi a a b
Laboratoire de Physique Théorique de la Matière Condensée (UMR CNRS 7600), Université Pierre et Marie Curie, 4 Place Jussieu, F75252, Paris cedex 05, France Department of Physics, Faculty of Sciences, University of Split, Nikole Tesle 12, 21000, Split, Croatia
a r t i c l e
i n f o
Article history: Received 3 February 2010 Received in revised form 20 April 2010 Accepted 7 May 2010 Available online 13 May 2010 Keywords: Aqueous mixtures Molecular liquids Simulations Statistical mechanics
a b s t r a c t Two binary mixtures with methanol as common component, namely water–methanol and methanol– acetone are studied by Molecular Dynamics simulations. Thermodynamical properties such as enthalpies, excess enthalpies, volumes and excess volumes are compared. Structural properties are studied through the various site–site pair correlation functions and the associated Kirkwood–Buff integrals. While the thermodynamical properties are relatively well calculated, we show that structural properties are affected by two severe problems. The first is an inherent property of the asymptote of the correlation function in finite systems, which affects their integrals, the second is the micro-heterogeneous structure of these hydrogen bonding mixtures, which affects the medium-to-long-range part of the distribution functions. These two problems conspire to make computer simulation unreliable to study this precise part of the structural properties of such mixtures. We propose pathways to correct this situation and demonstrate how it serves to considerably improve the calculated Kirkwood–Buff Integrals. Finally, we compare the two mixtures relative to their tendencies to form local heterogeneity. The analysis demonstrates that the microstructure of neat methanol is better preserved in acetone than in water. © 2010 Elsevier B.V. All rights reserved.
1. Introduction It is now relatively well documented and accepted that mixtures of associating liquids, such as aqueous mixtures, are difficult to simulate, and in particular, need excessively large simulation times and system sizes [1–5], which is a real paradox, mainly in view of the small size of these molecules. The reason for this unusual behaviour is that molecules preferentially self-associate through the strong hydrogen bonding interaction, and therefore, in addition to the molecular dynamical scale, they introduce a second dynamic scale — that of the clusters that are formed. As a result, one needs large system sizes, mainly to probe the statistics of the various clusters that are formed, and because of the slower motion of the clusters dynamics, one also needs excessively long run times. Paradoxically, this problem is less serious for the neat systems, such as water or alcohols, where, nevertheless, considerable micro-structure have been shown to exist [6–9]. The reason for this difference may be tracked back to the fact that clusters in neat systems have either a well defined geometry, such as those in neat alcohols [8], or are too elusive to be detected by ordinary statistics, such as water or amides [10]. In contrast, mixtures of such liquids tend to have some sort of self-partitioning [12,13], quite visible in snapshots, which confers these systems a rigidity that ⁎ Corresponding author. E-mail address: [email protected] (A. Perera). 0167-7322/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2010.05.006
defies simulation techniques. Nanostructuring is not unique to aqueous mixtures and can be observed in computer simulations of various ionic liquids as well [11], where hydrophobic chains selfsegregate from the ionic headgroups. There is an extensive computer simulation studies of aqueous methanol mixtures by several authors [3,12–17], and this system is now recognized to represent a typical example of well mixed system that exhibits strong micro-heterogeneity (MH). We revisit here this mixture, that we will consider as a reference system, in order to compare with the acetone–methanol mixture, which is supposed to be less structured, from the simple argument that methanol has only one donor and one acceptor bonding sites, as opposed to water which has two of each. In an earlier work [18], we have studied the acetone– water mixture and shown that most existing models lead to demixing schemes. In view of this unsatisfactory status, we study here the acetone–methanol mixture, which is well behaved with respect to the demixing issue. It is noteworthy that even this simple mixture is not satisfactorily simulated with usual force field generation techniques and specific issues have been discussed in recent works [15–17]. There are less studies of the acetone–methanol mixtures, to our knowledge. Computer simulation studies [19,20] of systems of less than 600 molecules report the existence of local heterogeneity and methanol chain clustering. This appears to be confirmed by various spectroscopic techniques [21–23], although the extent of segregation differs according to different techniques [22,23].
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
The main focus of this paper is how micro-heterogeneity influences the calculation of the structural functions and how to deal with this problem, and when this is possible to deal with it efficiently. For this reason, it is important to examine mixtures with expected different behaviour and see how they are affected by their respective micro-heterogeneity. 2. Simulation and theoretical details 2.1. Simulations In this work, we have considered OPLS force field models for both acetone [24] and methanol [25], and SPC/E model for water [26]. We emphasize again that neither SPC/E nor TIP4P water could mix satisfactorily with OPLS acetone [18]. The problem of detecting a true demixing within finite size computer simulation is not easily resolved when strong micro-heterogeneity is present. This important issue has not been satisfactorily addressed at present. The problematic is the following. Phase separation is induced by strong concentration fluctuations. These can be tackled properly within computer simulations by scaling techniques [27]. However, when local heterogeneity is present, it becomes difficult to disentangle this feature from genuine concentration fluctuations. For both mixtures, we have performed Molecular Dynamics simulation in the isobaric ensemble, for the ambient conditions of temperature and pressure. Berendsen thermostat and barostat where applied, with relaxation times of 0.1 ps and 0.5 ps, respectively. The time step for the integration of motion was fixed at 2.10− 15 s in all the calculations, which was found satisfactory in order to reproduce the rotational motion of water. In all our simulations total of N = 2048 molecules is considered for each mixture. This size allows us to reach intermolecular distance between 2.5 nm and 3 nm, which is about 8 to 10 times the molecular size, and which should be, in principle, sufficient to reach proper asymptotic bulk behaviour for the pair correlation functions. This is not often true, as we discuss in the results section. For each mixtures, the mole fraction of one of the component was varied by steps of 0.1, from zero to one. Each such system was equilibrated for 100–120 ps. This time window is sufficient to establish equilibrium values for the energies and volumes. Statistics were performed in batches of run length varying in the range 64– 256 ps, for about 10–20 batches. The accumulated statistics cover a range from 1 to 3 ns. This is still about twice smaller than what is claimed to be optimal in recent works [1]. The DLPOLY2 code [28] was used to perform these simulations. As pointed out in our earlier work [5], this code does not reproduce the expected asymptotic behaviour of the site–site correlation functions, and a multiplicative correction factor N/(N − 1) needs to be applied. This problem is independent of the issue of the asymptotical behaviour of these functions inherent to finite size, that we now expose. 2.2. Handling the asymptotic behaviour of pair correlation functions in mixtures Although computer simulations are performed with periodic boundary conditions, they are nevertheless not handling infinite size system. It is a well known result of statistical mechanics [29] that the asymptotical behaviour of the radial distribution function of an finite system, for example the canonical or isobaric ensemble, goes to a limit which is 1 − 1/N rather than exactly one as it should be for an open system (grand canonical ensemble). In the thermodynamic limit, which is when N (the number of particles) and V (the volume) go to infinity, then the equivalence of all these ensembles is automatically achieved [29]. But, in practice, the problem of the asymptote not being unity is a cumbersome one if one wants to compute the Kirkwood–Buff integrals (KBI), which are integrals over the pair distribution functions. Surprisingly, this important point has not been pointed out in many
53
earlier works concerning the evaluation of the KBI. But there is even more. In many textbooks, the limit that is proposed, which is 1 − 1/N, is in fact that of the ideal gas system. In presence of interactions and correlations, the proper result has been known for decades from the earlier works [30,31] of Lebowitz and Percus in particular, and it is 1 − 1/(NκT⁎), where κT⁎ = ( ∂ βP / ∂ ρ)T is the reduced isothermal compressibility, where β = 1/kBT is the Boltzmann factor, P the pressure and ρ the density. For the ideal gas, κT⁎ = 1, and the often quoted limit is recovered. Otherwise, there is a clear asymptote to below unity, that can create problems if the KBI are evaluated from this RDF. The extension of this expression for mixtures [32], in terms of the pair correlation between the pair of species i, j, reads: lim g ð1; 2Þ r→∞ ij
1 1 ∂ρ = 1− qffiffiffiffiffiffiffiffiffiffi ð i ÞTNk Ni Nj ρ ∂βμ j
ð1Þ
where μi is the chemical potential of species i. The evaluation of this derivative in practical calculations is a tedious task, since a chemical potential is needed, which is already a very noisy quantity to obtain from simulations, hence its derivative would be entitled to even larger errors. It is important to realize that, for the case of mixtures, the 1/N factor does not necessarily imply that the asymptotical value will be small if N is large. Indeed, for small concentration of one of the species, the corresponding factor NiNj can be arbitrarily small. We illustrate this in the results section. The conclusion of this fact, is that the asymptote of the pair distribution function can only be revealed by looking closely at the calculated functions. The expression above holds equally, as such, for the site–site correlation functions, since all site– site functions integrate out to give the same KBI as that of the centerof-mass correlation function [29]. This result is a consequence of the invariance of this quantity on the choice of the molecular center [29]. We propose a simple way to cure this problem artificially. We first spot the asymptotical behaviour of the function by direct visualization. Then, we add the following correcting function to the correlation function: " correct gij ðrÞ
=
raw gij ðrÞ
δij 1+ 2
r−κij 1 + tanh αij
!# ð2Þ
where the coefficients δij, κij and αij depend on the considered site–site correlation function. Usually, we take κij to be twice the distance at which gij(r) starts to be different from zero (an effective diameter). δij corresponds to the shift observed in the incorrect asymptote value aij with δij = 1 − aij, and αij controls the smoothness of the passage from aij to 1. We take this parameter to be a constant factor about twice the molecular diameter. We have verified that this value could be varied from the molecular diameter to about 3 times this distance without noticeable changes in the integrals. The reason of this invariance is due to the fact that the amplitude of the change is of the order of 0.1%. However small this value could be, it nevertheless leads to noticeable deviation of the KBI as will be illustrated in the next section. The form of the function above avoids perturbing the original function near the first few neighbours, and in particular in the core region, which is an absolute requirement. Finally, we note that, for neat systems, this tail correction can be so small, usually when system sizes larger than N = 1000 are considered, that the calculated RDF tend to unity within the numerical error [33]. This is particularly the case when the compressibility of the neat liquid is very small, as it is often the case with most liquids. 3. Results We present thermodynamic properties separately from the structural ones principally because the first ones do not seem to be affected the same way as structural ones. Properties such as enthalpies
54
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
and volumes converge much faster towards their equilibrium values and need only modest system sizes, such as N = 300–500 molecules. This is in sharp contrast with the proper computation of the radial distribution functions, which requires larger systems and long run times. This important fact indicates that the strong local inhomogeneity has a definitive impact on the evaluation of different properties. This is in contrast with the evaluation of the same properties in simple Lennard–Jones liquids and mixtures, where the statistics of both quantities require the same spatial and temporal scales, except near phase separation boundaries. 3.1. Thermodynamic properties We present the enthalpies and the volumes, as well as their excess quantities defined for the quantity A(x) as Aex(x)=A(x) −A1 −x(A2 − A1), as a function of x, the mole fraction of species labeled 2, and A1 =A(0) and A2 =A(1) are the neat liquid values of species 1 and 2, as obtained for x=0 and x=1, respectively. We compare our results with experiments and other works for the case of aqueous methanol solution. The reported enthalpies are in fact configurational energies, since we maintain atmospheric pressure, the PV term is negligible [34]. Fig. 1 shows the densities (top panel), volumes (middle panel) and excess volumes (lower panel) for the aqueous methanol mixture (labeled from now WM) while Fig. 2 shows the same quantities for the acetone–methanol mixtures (labeled AM). There is an obvious difference in the excess quantities: that of the AM mixture is nearly zero, indicating very little volume change when mixing the two components. Since mixing two hard spheres liquids leads to a noticeable negative excess volume [35], this means that the flat excess volumes of AM mixture indicate an increasingly repulsive excess volume as acetone is added into the mixture. Since hard spheres mix homogeneously, none of the tendencies observed here demonstrate any a priori peculiar behaviour related to the formation of micro-heterogeneity. We recall that
Fig. 1. The densities (top panel), volumes (middle panel) and excess volumes (lower panel) for aqueous methanol , as functions of the methanol mole fraction. The red curve is experimental data from Ref. [36] and excess volumes from Ref. [37]. Black symbols are for the present MD data, green symbols from simulations of Ref [3] and magenta symbols from Ref. [14]. In the lower panel blue symbols are from Ref. [17].
Fig. 2. The densities (top panel), volumes (middle panel) and excess volumes (lower panel) for acetone–methanol mixture, as functions of the acetone mole fraction.
the aqueous acetone excess volumes were found negative for many model combinations while most showed tendency for demixing [18]. Fig. 3 shows the energies and excess energies for the WM mixture and Fig. 4 for the AM mixture. The principal difference between the two mixtures is that the excess for WM are negative while they are positive for the AM mixture. This may be thought to be consistent with the behaviour found for the excess volumes: the reason for near flat excess
Fig. 3. The energies (top panel) and excess energies (lower panel) for aqueous methanol, as functions of the methanol mole fraction. Line and symbols as in Fig. 1. The red line is experimental data from Ref. [38].
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
55
MH problem, and then by computing the KBI, which are thermodynamic quantities that relate directly to probing through measurable quantities such as partial molar volumes and chemical potentials. 3.2.1. Snapshots Fig. 5 shows a snapshot of the aqueous methanol mixture for 30% methanol content, where the water(left) and methanol(right) molecules are shown separately for the same angle of shot. It is visually obvious that water molecules have a Swiss-cheese like distribution, noticeable from the fact that distant molecules are seen as smaller than those closer to the viewer. This is particularly obvious when the shots are made to rotate when visualizing them. In contrast, the methanol molecules appear more homogeneously distributed. This means that some of the water clusters contain embedded methanol molecules. Similar snapshot-based conclusions about micro-segregation were reported in other works [12] and for other mixtures [1,10,18]. Fig. 6 shows snapshots of the acetone–methanol mixture for 30% acetone content. The micro-heterogeneous distribution of the methanol molecules (left) is quite noticeable, in contrast to that of the acetone molecules (right) which appears more homogeneous. These observations indicate that the common species, methanol, plays a different role in different environment. In particular, these shots seem to imply that it is always the more strongly associating molecule that drives the micro-segregation. Fig. 4. The energies (top panel) and excess energies (lower panel) for the acetone– methanol mixture, as functions of the acetone mole fraction.
volumes of the AM mixture can come from the overall repulsiveness witnessed through the excess energies. We recall that the excess energies for the aqueous acetone mixtures were found positive for nearly all mixtures [18]. In contrast, however, the AM mixture appears well mixed, despite the inherent micro-heterogeneity. These various trends do not allow to guess the micro-heterogeneous status of these mixtures, similarly to that found for the volumes. We note that, for the aqueous methanol mixture, most calculated energies are in relatively good agreement with the experimental ones, except for the fact that all calculated excess energies have a minimum at x = 0.5 while experimental ones are positioned at x ≈ 0.3. This strange fact has not been yet satisfactorily interpreted. 3.2. Structural properties The study of these properties can be made in three steps. First by showing evidence of MH through the snapshots, next by studying the correlation functions themselves, which is a statistical approach to the
3.2.2. Distribution functions Fig. 7 shows selected site–site radial distribution functions (ssRDF) to illustrate the structure of the aqueous methanol mixture for methanol mole fractions of 0.20, 0.50 and 0.80. These plots indicate that the sensitivity of the water–water and meth–meth correlations, to the variation of methanol mole fraction, is quite different. Water– water correlations (in blue) show a dramatic increase of the first few neighbours structure, which indicates that water forms tight clusters in methanol rich environment. In comparison, the various methanol– methanol correlations seem relatively insensitive to changes in water concentration. In particular, the methanol oxygen–oxygen (O–O) correlations show the inverse trend to that of water: they increase with methanol concentration while the water–water correlations increase with decrease of water concentration. Since both species are hydrogen bonding species, one would expect a symmetric behaviour, and the observed asymmetry is rather striking. The interpretation of this fact corroborates the snapshots of Fig. 5. Water tends to microsegregate more than methanol, resulting in an increased correlation when it is scarce in methanol, as it tends to form small tight clusters. Conversely, methanol is more homogeneously distributed in water, as
Fig. 5. Snapshots of the aqueous methanol mixture at 30% methanol content. Water molecules are shown in the left and methanol molecules at the right. In both shots the missing molecules are shown in skeleton. The color code is white for H atoms, red for O atoms and blue for the methyl group.
56
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
Fig. 6. Snapshot of the acetone–methanol mixture at 30% acetone content. Methanol molecules on the left and acetone molecules on the right. The color codes are the same as in Fig. 5.
can be seen from the little variation of the methyl group correlations. There must be some micro-segregation of methanol at low content (20%), since its distribution is nearly that of 80% methanol content, but one cannot really speak of methanol cluster the same way water is clustered in 80% methanol, because of the marked differences between the first peaks in the respective cases. By comparing the structure between the different meth–meth correlations we can guess the way methanol molecules orient in various conditions. In neat methanol, we have previously compiled the structure of different methanol clusters [9,32] and noted that the methyl (M) sites were homogeneously distributed while the O sites were connected in
Fig. 7. Selected site–site distribution functions for the aqueous methanol mixture, at methanol mole fraction x = 0.20, 0.50 and 0.80. The oxygen–oxygen (OO) are shown in blue for water, black for methanol and magenta for water–methanol. The red curve is for methanol MM correlations and the green curve for water–Methanol OO cross correlations.
various chain-like patterns. Assuming this is still the case in 80% meth content, we see that the decrease of the O–O correlations with methanol content indicates that methanol O sites would preferentially bond with water. The increasing cross correlations between the oxygens of the two species (in magenta), with methanol decrease, confirm this hypothesis. In contrast, the relative insensitivity of the M–M correlations (in red) indicate that the methyl sites are homogeneously distributed within the methanol rich areas in all circumstances. These various trends confirm the local micro-heterogeneous state of the mixture. Fig. 8 shows selected ssRDF for the acetone–methanol mixture. Here, only methanol is self-bonding. Therefore, it is not surprising to see methanol O–O correlations (in magenta) behave similarly to water–water correlations in Fig. 7. Indeed, the first neighbour correlation increases dramatically with acetone concentration, due to strong clusterisation of methanol in acetone rich environment. In comparison, we see that the acetone O–O correlations (in blue) are insensitive to changes in concentration, indicating that acetone is correlated to itself much the same way in concentrated and diluted situations. This is in favour of micro-segregation of the two species. The methanol M–M correlations (in red) increase when methanol decreases, although less dramatically than the O–O correlations. This indicates that the M sites are not so randomly distributed (as they appeared to be in aqueous methanol), enforcing the picture that the methanol clusters should be quite similar to that of the neat methanol (chain-like). The O–O cross correlations (in green) indicate that little Hbonding occurs between the two species. These two figures illustrate the difference in micro-heterogeneity in the two types of mixtures sharing a common Hbonding molecule. From the detailed analysis above, the major difference is in the cluster structure of methanol in the two different environments. The microstructure of neat methanol clusters appears to be more preserved in acetone than in water. We turn now to the asymptotical behaviour of these distribution functions. Fig. 9 serves to illustrate the fact that the asymptote of the correlation functions do not go to unity, but behave rather like in Eq. (1). For this purpose we pick the acetone–methanol mixture at 20% acetone content, and plot the following two methanol–methanol ssRDF, oxygen–oxygen and methyl–methyl. The original functions are shown in dashed lines and those corrected for the asymptote in full lines. The top panel shows these 2 functions and the differences with the corrected ones cannot be distinguished on this scale. The middle panel shows a zoom on the tail part of the ssRDF. The curve in red shows the correcting function suggested in Eq. (2), and the green curve shows the current asymptote, which we estimate at a = 0.9983 from the plot, which is obviously different from unity. It is important to note that both ssRDF tend to the same asymptote. The curves
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
Fig. 8. Selected site–site distribution functions for the acetone–methanol mixture, at acetone mole fraction x = 0.20, 0.50 and 0.80. The magenta curve is for meth–meth O–O correlations, the red curve for meth–meth M–M correlations, the blue curve for acetone–acetone O–O correlations and the green curve for acetone–methanol O–O cross correlations.
corrected for the asymptote are shown in full lines, with the following parameters in Eq. (2), κ = 5Å and α = 1 Å. The lower panel shows the corresponding running KBIs (RKBI) defined as [10] r
2
Gij ðrÞ = 4π∫0 dt t ðgij ðrÞ−1Þ
ð3Þ
where gij(r) is the ssRDF between sites i and j. The asymptotic behaviour of this latter function, when r → ∞ is precisely the KBI. There is a major difference between corrected and uncorrected KBIs, so the correction is really needed. The downwards curvature of the uncorrected curve (in dashes) comes simply from the integrated r3 behaviour in Eq. (4), that arises if the asymptote is not strictly unity. It is important to note that this correction is not due to the existence of MH, but to the finite size effect, and can also be seen in simple liquids [31] and their mixtures [32]. However, the existence of MH can mask this behaviour, depending on the extension of correlations which are affected by the fact that this micro-heterogeneity extends up to the half box size, which is the usual cutoff-distance for the computation of the correlations. Finally, we show an example of the behaviour of the correlation function where the existence of MH does not allow a proper settling of the asymptotical behaviour of the RDF and hence of the corresponding RBKIs. This is illustrated in Fig. 10, where we show selected ssRDF for the equimolar aqueous methanol mixture. Each curve presents the outcome of a batch of statistics for run length of 64 ps. The average that we consider acceptable is shown in thicker black line. It is clear that there is no convergence of the behaviour to any proper definite value. This type of problem is recurrent in the simulation of any aqueous mixture, and is never reported in the literature. It should also occur for other microsegregating functions, such as for aqueous ionic-liquid mixture for
57
Fig. 9. (Top panel) The oxygen–oxygen (OO in blue) and methyl–methyl group (MM in magenta) site–site correlation functions for the acetone–methanol mixture with 20% acetone. (Middle panel) Zoom on the asymptotical behaviour of the OO and MM correlation of the top panel. The correcting function (Eq. (2)) is shown in red (parameters in the text). The current asymptote is shown in green. The corrected functions are shown in full lines and the uncorrected ones in dashes. (Lower panel) The RKBIs corresponding to the functions above, with the same color conventions.
example. These small differences that are amplified in the evaluation of the RKBI are unnoticeable in the first few peaks of the full RDF (upper panel). This could be the reason why these problems are never noticed in the published literature, when showing the RDF under the 12 Å range. We wish to emphasize that correlations at the first and second peak are relatively well stabilized after few 100 ps, and it is only the behaviour beyond that is subject to the anomalous behaviour reported here. Hence, results reported for smaller system sizes should give a reliable account of the correlations of the first two neighbour shells. But one should keep in mind that this fact masks the inherent pathological behaviour for higher neighbour shells, that is ultimately due to a real physical effect, that of the micro-heterogeneity. 3.2.3. Kirkwood–Buff integrals In Fig. 11, we report the KBI of the two mixtures, together with the experimental results for the aqueous methanol mixture. We are not aware of any experimental results for the acet–meth mixture. In place, we have plotted an estimate of the experimental KBIs, evaluated from the following expressions [32]: G12 ðxÞ = −
V1 V2 VðxÞDðxÞ
Gii ðxÞ = G12 ðxÞ +
Vj 1 −VðxÞ xi DðxÞ
ð4Þ
where x = x2 is the mole fraction of component 2 (methanol in water– meth mixture and acetone in acet–meth mixture), V1 and V2 are the pure component volumes, V(x) = V1 + x(V2 − V1) is the approximate
58
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
Fig. 10. (Top panel) The water–water radial distribution function of aqueous methanol equimolar mixture for different statistics (undistinguishable on this scale). (Middle panel) Zoom on the asymptotical behaviour of the RDFs above, each curve represents a statistics over 64 ps. (Lower panel) The RKBIs corresponding to the functions above, with the total average shown in the thicker black curve.
volume of the mixture and D(x) = 1 − x(1 − x)d is related to the chemical potential derivative in a crude way, and i = 1,2 with j = 2,1 in the second formula. These expressions, which neglect the compressibility of the mixture as well as specificity of the partial molar volumes and chemical potentials, reproduce quite well the experimental KBIs of the aqueous methanol mixture. It is often found that modifying the definition of D(x) can allow a very good fit of the KBIs [39]. The volume for the neat water is taken as V = 18 cm3/mol, for methanol V = 40 cm3/mol and for acetone V = 78 cm3/mol. We use d = 1.20 for aqueous methanol and d = 2.10 for the acetone–methanol mixture. The MD results are obtained by correcting the asymptotic result whenever this is possible. As noted above, this is possible when the asymptote is not clouded by bad statistics due to MH. We find that this was invariably the case when considering self-correlations of the majority species, as illustrated in Fig. 9. Conversely, self-correlation of the minority species was always difficult to correct since the statistics of the clustering strongly affect the stabilisation of the asymptote. In such bad cases, the asymptotical behaviour is chosen by averaging the different values (see Fig. 10). The reported results indicate that the KBI evaluated after correction are nevertheless in very good agreement with the experimental results for the aqueous methanol mixture and with the fitting expressions of Eq. (4) above. The first thing that we notice is that the KBIs are relatively well defined for the solute (component 2) in almost the full solute mole fractions, while for the solvent this is true only in the solvent rich region. In general, the component poor region leads to poor estimate of the corresponding KBIs. This is clearly due to the large clustering of such species in the adverse rich component, which in turn affects the behaviour of the tail of the ssRDF. We note that methanol in the two mixtures has a similar behaviour when it is in low concentration: the KBI are noisy and out of bounds. This is again due to clustering and therefore bad definition of the asymptote. We note again that methanol, that is the common species, behaves very differently when it is in water than when it is in acetone. The striking similarity of the behaviour of methanol in aqueous methanol with that of acetone in acetone–methanol, as well as that of methanol in acetone–methanol with that of water in aqueous methanol, emphasizes the asymmetrical role of this substance in different environments. These particularities are due to the Hbonding nature of the species. The next obvious remark is that there is no direct signature of the differences in the micro-heterogeneous behaviour of these two mixtures that can be detected from simple visualization of the KBIs. Therefore, the ssRDF remain the best observables to appreciate the MH of a mixture, albeit the problems encountered in evaluating them. 4. Discussion
Fig. 11. The KBI for the aqueous–methanol mixture versus methanol mole fraction (top panel) , and acetone–methanol mixtures versus acetone mole fraction(lower panel). MD results are shown in filled dots (top panel: blue for water–water, red for Meth– Meth and black for Water–Meth; bottom panel: blue for Meth–Meth, red for Acet–Acet and black for Meth–Acet). Experimental results from Ref. [39] are shown in the top panel in open squares with color convention similar to the MD symbols. The fit from Eq. (4) are shown in lines with color convention similar to the MD symbols.
The fact that the simple asymptote correction gives appreciably good KBI indicates the importance of such correction. It is also obvious that rather large system sizes must be studied, as can be clearly seen from the zoom of Figs. 9 and 10: one needs to reach the part where oscillations due to packing structure die out, and the bulk regime starts to show up. Generally, this happens roughly after 5 molecular diameters, although no general rules can be easily specified. The paradox we pointed out earlier, namely that one needs such large systems and therefore long statistics and runs for such small molecules mixtures, is truly disconcerting one. We would like to stress that these calculations were conducted with unmodified versions of the OPLS force fields of the methanol and acetone molecules. The fact that a good expectation of the KBIs is obtained in this work, enforces the idea that it may not be really necessary to modify existing force fields, such as the OPLS force field for example, in order to achieve this goal. Such alternate routes were indeed taken recently with success by some authors [1–4]. The present study
A. Perera et al. / Journal of Molecular Liquids 159 (2011) 52–59
shows that the painstaking task of redesigning force fields could be avoided by explicitly taking into account inherent corrections due to finite size effects. Following this route, it seems that such type of corrections could be extended to account for the spurious contributions due to the local heterogeneity. We would like to stress, however, that in cases where spurious demixing is enforced by improper force fields [18], it would be rather difficult to use such tricks to improve results on KBI, for obvious physical reasons. Ultimately, the pathways followed in this work amount to distinguish between contributions due to concentration fluctuations and those due to local micro-heterogeneity [5]. As stated in the Introduction, the first contributions control phase separation while the second are responsible for the difficulties encountered in simulations of these particular systems. It is at present quite difficult to predict how these two contributions stem from the force field parametrisation. Recently, some authors [16,17] have raised the subtle issue of reproducing the experimentally observed minimum in the excess partial molar volume of the aqueous alcohol mixtures. It seems that none of existing force fields are able to reproduce such minimum. Another, perhaps related issue, is the shift in minimum in the experimental excess enthalpy of aqueous methanol around x = 0.28, instead of 0.5 for all force field models (see lower panel of Fig. 3). Ultimately, a convincing computer simulation study would have to be able to account for the subtle local structural orders that can give rise to these still unexplained features. 5. Conclusion This paper focuses on the micro-heterogeneous structure of different types of hydrogen bonding mixtures and the ability to report on this property through MD simulations. For this purpose, two mixtures have been studied with a common component, methanol. The main conclusion of this work is that these micro-segregating mixtures can be studied by MD simulation in a reliable way if one focuses on properties such as enthalpies and volumes, as well as correlation functions when these are limited to the first and second neighbours shells. Generally speaking, these correlations tend to die rather quickly after the second and third neighbours range. This is why most reported results in the literature are never concerned with this problematic. However, this apparent success should not mask the real underlying problem caused by the fact that these mixtures are strongly inhomogeneous, and that such inhomogeneity affects undeniably long range properties of these mixtures in an uncontrollable way. This is due to two facts: first the inherent problem of the asymptote, for which we provide here a remedy, and second the extent of the local clustering that can affect the convergence of the value of the asymptote. Increasing the size of the system may appear as a remedy, but it comes with additional
59
problems such as the need for longer runs, since one should now also take into account the slow dynamics of the clustered parts of the simulated sample. The present calculations show that methanol behaves differently in the aqueous and acetone environments, mainly with respect to clustering. It clusters less in water than in acetone. This fact alone indicates that micro-heterogeneity is not an intuitively predictable property of these mixtures.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]
R. Chitra, P.E. Smith, J. Phys. Chem. B104 (2000) 5854. R. Chitra, P.E. Smith, J. Chem. Phys. 115 (2001) 5521. S. Weerasinge, P.E. Smith, J. Chem. Phys. 118 (2003) 10663. M.E. Lee, N.F. van der Vegt, J. Chem. Phys. 122 (2005) 114509. L. Zoranic, R. Mazighi, F. Sokolic, A. Perera, J. Chem. Phys. 130 (2009) 124315. J.-H. Guo, Y. Luo, A. Augustsson, S. Kashtanov, J.-E. Rubensson, D.K. Shuh, H. Agren, J. Nordgren, Phys. Rev. Lett. 91 (2003) 157401. R. Ludwig, ChemPhysChem 6 (2005) 1369. L. Zoranic, F. Sokolic, A. Perera, J. Chem. Phys. 127 (2007) 024502. A. Perera, F. Sokolic, L. Zoranic, Phys. Rev. E75 (2007) 060502-(R). L. Zoranic, R. Mazighi, F. Sokolic, A. Perera, J. Phys. Chem. C111 (2007) 15586. J.N.A. Canongia Lopes, A.A.H. Padua, J. Phys. Chem. B110 (2006) 3330. S. Dixit, J. Crain, W.C. Poon, J.L. Finney, A.K. Soper, Nature 416 (2002) 829. L. Dougan, S.P. Bates, R. Hargreaves, J.P. Fox, J. Crain, J.L. Finney, V. Reat, A.K. Soper, J. Chem. Phys. 121 (2004) 6456. A.M. Ferrario, M.M. Haughney, I.R. McDonal, M.L. Klein, J. Chem. Phys. 93 (1990) 5156. E.J.W. Wensik, A.C. Hoffmann, P.J. van Maaren, D. van der Spoel, J. Chem. Phys. 119 (2005) 7308 52003. L. Vlček, I. Nezbeda, J. Mol. Liq. 131–132 (2007) 158. Diego González-Salgado, I. Nezbeda, Fluid Phase Equilib. 240 (2006) 161. A. Perera, F. Sokolic, J. Chem. Phys. 121 (22) (2004) 11272. D.S. Venables, C.A. Schmuttenmaer, J. Chem. Phys. 113 (2000) 3249. G. Kamath, G. Georgiev, J.J. Potoff, J. Phys. Chem. B109 (2005) 19463. D.S. Venables, C.A. Schmuttenmaer, J. Chem. Phys. 113 (2000) 3243. M. Musso, M.G. Giorgini, H. Torii, J. Mol. Liq. 147 (2009) 37. J.-J. Max, C. Chapados, J. Chem. Phys. 122 (2005). W.L. Jorgensen, J.M. Briggs, M.L. Contreras, J. Phys. Chem. 94 (1990) 1683. W.L. Jorgensen, J. Phys. Chem. 90 (1986) 1276. J.C. Berendsen, J.P.M. Postma, W.F. Von Gusteren, J. Hermans, in: B. Pullman (Ed.), Intermolecular Forces, Reidel, Dortrecht, 1981. N.G. Wilding, Phys. Rev. E. 52 (1995) 602. DL_POLY package from CCLRC Daresbury Laboratory www.ccp5.ac.uk/DL_POLY/. J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986. J.L. Lebowitz, J.K. Percus, Phys. Rev. 122 (1961) 1675. J. Kolafa, S. Labik, A. Malijevsky, Mol. Phys. 100 (2002) 2629. A. Perera, L. Zoranic, F. Sokolic and R. Mazighi (in preparation). L. Zoranic, F. Sokolic, A. Perera, J. Mol. Liq. 136 (2007) 199. F. Sokolic, A. Idrissi, A. Perera, J. Chem. Phys. 116 (2002) 1636. J.S. Rowlinson, F.L. Swinton, Liquids and Liquid Mixtures, Butterworths, 1982. Landolt-Börnstein, Neue Serie Group IV , vol.1 part B, SPringer Verlag, 1977. Y. Koga, Can. J. Chem. 77 (1999) 2039. L. Benjamin, G.C. Benson, J. Am. Chem. Soc. 67 (1963) 858. A. Perera, F. Sokolic, L. Almasy, Y. Koga, J. Chem. Phys. 124 (2006) 124515.