On the composition dependence of thermodynamic, dynamic and dielectric properties of water-methanol model mixtures. Molecular dynamics simulation results

We have investigated thermodynamic and dynamic properties as well as the dielectric constant of water-metha\-nol model mixtures in the entire range of composition by using constant pressure molecular dynamics simulations at ambient conditions. The SPC/E and TIP4P/Ew water models are used in combination with the OPLS united atom modelling for methanol. Changes of the average number of hydrogen bonds between particles of different species and of the fractions of differently bonded molecules are put in correspondence with the behavior of excess mixing volume and enthalpy, of self-diffusion coefficients and rotational relaxation times. From the detailed analyses of the results obtained in this work, we conclude that an improvement of the description of an ample set of properties of water-methanol mixtures can possibly be reached, if a more sophisticated, carefully parameterized, e.g., all atom, model for methanol is used. Moreover, exploration of parametrization of the methanol force field, with simultaneous application of different combination rules for methanol-water cross interactions, is required.


Introduction
In the very recent work from this laboratory we reported a detailed analysis of molecular dynamics computer simulation data for the microscopic structure of water-methanol mixtures in the entire composition range, and performed comparisons with experimental results obtained using X-ray diffraction [1]. Our principal focus in that publication was on the changes, with composition, of the experimental and theoretical total structure factors of mixtures considered at room temperature and ambient pressure. On the molecular dynamics (MD) side, we analyzed the SPC/E [2] and the TIP4P/Ew [3] water models combined with the OPLS/AA (all atom) methanol model [4]. The all atom modelling for methanol was used in that study because it is intrinsically required by the procedure to perform comparisons with the experimental total structure factors.
Encouraged by an overall satisfactory performance of the MD simulation data in the description of the experimental trends, in the present work our primary objective is to extend our previous study by exploring the results for a more comprehensive set of properties of water-methanol mixtures, besides the microscopic structure. Namely, we focus on thermodynamic, dynamic and dielectric properties obtained from our MD simulations and perform comparisons with the available experimental data. Results of computer simulations from other authors are included as well, aiming at a more comprehensive insight into the dependence of properties of mixtures in question on the conditions of computer modelling.
There have been very many studies of water-alcohol mixtures, and of water-methanol mixtures in particular, with the use of computer simulation methods since the pioneering works by Jorgensen and Madura [5] and by Okazaki et al. [6,7], principally focused on "infinite dilution" solutions by considering a single methanol molecule in a set of water molecules. The results were obtained using Monte Carlo (MC) simulations in the NPT and NVT ensemble, respectively. To our best knowledge, Ferrario et al. [8] and Tanaka with Gubbins [9] were the first to report thermodynamic and dynamic properties of these systems in the entire composition range from MD and MC simulations, respectively. Since then, much knowledge has been accumulated about the properties of water-methanol mixtures, see e.g., [10][11][12][13]. In particular, one of the more recent works on the topic, with the principal focus on dynamic properties of wateralcohol mixtures (methanol, ethanol and 1-propanol) using the TIP4P water model [14] in combination with the OPLS/AA model for alcohols [4] has been published by Wensink et al. [15]. On the other hand, Guevara-Carrion et al. [16] explored the SPC/E and TIP4P-2005 [17] water models in combination with a united atom type methanol model [18] for mixtures, reporting excess mixing volume and enthalpy, self-diffusion coefficients, shear viscosity and power spectra. Recent works on water-methanol mixtures include detailed studies of the mixing volume and mixing enthalpy [19,20], using the NPT MC simulations.
Analyses of various results emerging from computer simulations by using non-polarizable water models and different versions of force fields for alcohols in comparison with experimental data clearly show the necessity of further work. Namely, it seems necessary to reach a systematic, more profound understanding of strengths and drawbacks of different models and provide a design aiming at a better agreement with experiments, desirably for an ample set of properties of mixtures of interest, rather than for a single particular property.
In this respect, in the present work we would like to report thermodynamic properties of mixing in terms of excess mixing volume and excess mixing enthalpy, a set of dynamic properties that include the self-diffusion coefficients and rotational relaxation times, as well as the dielectric constant, all together with the detailed analysis of hydrogen bonding within species and between them, in the entire range of composition, starting from pure water and terminating by pure methanol. Whenever possible, we discuss the results of other authors and make comparisons with experimental data from literature. We restrict our attention to two different water models and the OPLS/UA (united atom) potential model of methanol in order to establish their performance. The OPLS/UA model for methanol and other more complex alcohols is popular in comparison with the more detailed OPLS/AA model due to computational efficiency.
Several recent efforts have focused on the parametrization of different properties of water, see e.g., [21,22], targeting to improve the description of the dielectric constant and of the density anomaly. A similar approach has been applied to methanol at the united atom level, see [23], targeting the dielectric constant in particular. The strategy of the proposed parametrization includes scaling of site charges with the dielectric constant as a target, and scaling of the energy of non-bonded interactions leading to a corrected surface tension as subsequent steps. The critical temperature is used for control as well. The parametrization can be done in cycles to yield correct values for the targets. On the other hand, Schnabel et al. [18] developed a successful parametrization of the OPLS/UA type model for methanol using vaporliquid equilibrium data, focusing on the correct description of the hydrogen bonds statistics to correctly reproduce the nuclear magnetic resonance spectroscopic data.

Simulation details
All our simulations of water-methanol mixtures models were performed in the isothermal-isobaric ensemble at ambient pressure. The present work involves exploration of two commonly used water models, namely the SPC/E [2] and the TIP4P/Ew ones [3]. On the other hand, the OPLS/UA (united atom) rigid nonpolarizable model for methanol [24] in the present work is taken with the adjusted parameters of [25]. The OPLS/UA model for methanol has been designed so that the parameters of the nonbonded potentials between different types of atoms satisfy the geometric combination rules [24]. However, in the work of van Leeuwen and Smith [25], new parameters for nonbonded potentials were designed and standard

43602-2
Lorentz-Berthelot (LB) combination rules, i.e., ǫ i j = (ǫ ii ǫ j j ) 1/2 and σ i j = (σ ii + σ j j )/2, were used to calculate the liquid-vapor coexistence curve of methanol and a bit later of other alkanols [26]. Again, the LB combination rules were used for methanol in [27]. This is a setup similar to the one recently used by Moučka and Nezbeda [20], and we just apply the SPC/E and TIP4P/Ew models instead of the TIP4P potential considered by these authors. Moreover, we have not done any deviation from the LB combination rules, in contrast to [20]. All the simulations of water-OPLS/UA version of the models, unless specified, were performed using the DLPOLY Classic package [28]. We used the Berendsen thermostat and barostat with τ T = 0.1 ps and τ P = 2.0 ps, the running timestep was 0.002 ps. As common, periodic boundary conditions were used. The nonbonded interactions were cut-off at 11 Å, whereas the long-range interactions were handled by the Ewald method with a precision of 10 −5 . In order to maintain the geometry of water and methanol molecules, the SHAKE algorithm was used. For all simulations, the initial configurations were constructed by using from 2000 to 3000 molecules of the two species placed randomly in a cubic box. After equilibration, several sets of simulation runs were performed, each lasting for 6-8 ns, with restart from the previous configuration, to obtain the averages for data analysis. Actually, the overall length of simulation was dictated by the stability of internal energy, density and the dielectric constant, so that the runs were not shorter than 20 ns. The precise numbers of molecules of both species used in our simulations are given in table 1. Changes of the mixture compositions (as well as the compositions themselves) throughout this study are described by X m (molar fraction of methanol); its values are given in table 1 as well. In addition, we provide the weight concentration of methanol for convenience of the reader.

Results and discussion
First we would like to discuss our observations concerning the properties that do not require any reference to the species forming the mixture (such as overall density, mixing properties and the dielectric constant) and next proceed to the properties characterizing individual species characteristics such as self-diffusion coefficients and reorientational relaxation times, and discuss hydrogen bonding in terms of the average number of bonds per molecule of the particular species, cross bonds and bonding states of molecules of each species.

Density, mixing properties and the dielectric constant
Changes of the liquid mixture density with methanol mole fraction for the two combinations of potential models are shown in figure 1. The NPT molecular dynamics simulation data at room temperature (298.15 K) and at ambient pressure are accompanied by the experimental values [29]. The SPC/E-OPLS/UA model very well reproduces the experimental results. Only at high methanol fractions, the mixture density is slightly overestimated in comparison with the experimental data but this disagreement is hardly visible on the scale of the figure. It is worth mentioning that if the methanol component is described within the OPLS/AA model, then the density of the mixtures in question is described less satisfactorily: it is underestimated practically in the entire composition range [1]. To conclude with this figure, it seems that the united atom model for methanol is reliable to describe changes of the density on composition of the mixture at room temperature and at ambient pressure. However, a more sensitive test for the species modelling is illustrated by the calculations of the excess mixing volume, ∆V mix , figure 2. It is obtained straightforwardly as: the united atom level of modelling, changes of the water model yield a marginal improvement of the composition dependence of ∆V with respect to the experimental data [30]. Actually, the results are reasonably good for both water-rich and methanol-rich mixtures compared to the experiment. However, at the intermediate region, the applied modelling underestimates the nonideality of mixtures. Similar trends (to underestimate nonideality) have been observed for the TIP4P-OPLS/UA combination of models, if the geometric mixing rules were used for the cross-energies and cross-diameters [19], cf. figure 1 (a) of reference [19]. Moreover, the overall performance of the excess mixing volume coming from simulations does not improve for different explored sets of deviation parameters from the standard LB combination rules as shown by [20] for the TIP4P-OPLS/UA model, cf. figure 2 of reference [20]. Our figure 2 complements the results shown in figure 15 of reference [16]. Summarizing our results and insights coming from the previous simulations in [16,19,20], it seems that different combinations of water and methanol models yield a satisfactory description of the excess mixing volume and correctly predict the composition value, X m ≈ 0.5 at which ∆V mix is minimum. However, it is doubtful whether such a moderate success would be preserved for other temperatures and pressures and/or for other water-alcohol mixtures. Therefore, accumulation of knowledge about the performance of different models under other conditions and for other members of water-alcohol mixtures family is of great importance. On the other hand, it seems nec-

43602-4
Composition dependence of the properties of water-methanol mixtures essary to explore parametrization and/or sophistication of the OPLS/AA modelling, placing the focus on this particular experimentally measurable property.
The excess enthalpy of mixing, ∆H mix is determined as follows: where H mix is U pot + PV . As we see from the curves in figure 3 and from experimental results [31], the behavior predicted by simulations with either water models of this study, combined with OPLS/UA, is not satisfactory. On the methanol-rich side, the behavior of the theoretical curves seems to be a little better than on the water-rich side. However, in the entire range of compositions, computer modelling yields much weaker effects of nonideality compared to the experimental data. Moreover, the minimum of ∆H mix from simulations is seen at a higher methanol fraction compared to the experimental results.

43602-5
ence [19]. This behavior points out that the energetic aspects of mixing methanol and water in simulations should be reconsidered. A growth of the absolute value for ∆H mix can be reached by employing certain deviations to the combination rules [20], but then the position of the minimum for ∆H mix shifts to methanol-rich mixtures, which is in contradiction to the experimental behavior. Thus, it does not seem that the combination rules are responsible for this type of failure. One can hope that modelling of methanol at the all atom level would be profitable to improve an agreement with experimental data for this particular property. However, another possibility is that nonpolarizable models are intrinsically incapable of providing a desired agreement for this particular property. Definitely, other characteristics of the system that involve potential energy, besides excess mixing enthalpy, should be explored to make stronger conclusions. Now we proceed to another property that, as seen below, is not very successfully described within the framework of the potential models applied. It is known that the long-range, asymptotic behavior of correlation functions between molecules possessing dipole moments is described by the dielectric constant, ε. In general terms, the calculation of the dielectric constant from simulations is a demanding task, see e.g., [21,32,33] for the discussion of calculations of ε for pure water. Several factors can influence the result, such as the number of molecules, type of thermostat and barostat, precision of the summation of long-range interactions. Moreover, long runs are necessary to obtain reasonable estimates for this property, since it is usually calculated from the time-average of the fluctuations of the total dipole moment of the system [34] as follows: where k B is the Boltzmann constant and V is the simulation cell volume.
Both of the models for the pure species, water and methanol, are perfect in predicting their respective dielectric constants. As a result, the behavior of mixtures in the entire composition range is not perfect either, figure 4 (a). The inclination of the curve for ε(X m ) is not bad at all, but the values are too low compared to the experimental results. We get confidence in the quality of the present results making long simulations [ figure 4 (b)] with a sufficiently large number of particles in the simulation box. A simultaneous improvement of water and methanol models is necessary to obtain better values of the dielectric constant for the mixtures in question. How good is the parametrization of the methanol model of [18] with respect to the dielectric constant is unknown so far. On the other hand, quite recently the united atom model of methanol was parameterized to reproduce the experimental value of the dielectric con-

43602-6
stant [23], at the expense of making other properties worse, e.g., the self-diffusion coefficient. However, the performance of this modified model for water-methanol mixtures has not been evaluated yet. We plan to improve the description of the dielectric constant by performing simulations for more sophisticated models. Without such efforts it is difficult to attempt a successful description of the thermodynamics and the dynamics of ionic solvation in mixed water-alcohol solvents.

Dynamic properties and hydrogen bonding
We proceed to the description of our results for a set of dynamic properties. The self-diffusion coefficients for the two species of the mixture were calculated by two routes. Namely, from the mean-square displacement (MSD) of a particle via the Einstein relation, where i refers to water or methanol species. The MSD procedure to obtain D i was applied by using particle coordinates coming from DLPOLY simulations. On the other hand, we also obtained the selfdiffusion coefficients from the time integral of the velocity autocorrelation functions, For this purpose we employed the GROMACS software for the same model and with the same technical setup as in the DLPOLY package. This has been done due to the speed efficiency of GROMACS. For the sake of checking mutual consistency, a set of calculations has been performed and the results are given in figure 5. We observe that the two procedures yield very similar results, margins of inaccuracy overlap for all compositions. The subscripts "m" and "w" refer to methanol and water, respectively. A set of results for self-diffusion coefficients coming from simulations by using VACF is given in figure 6. While an overall behavior of the curves for water and methanol molecules is correct, the experimental [36] and simulation data differ in several details. On the water-rich side, the simulation data overestimate the self-diffusion coefficients as a result of inability of the explored models for water to yield the self-diffusion coefficient in perfect agreement with the experimental value for pure water. The self-diffusion coefficient for pure methanol is not as bad. The OPLS/UA model slightly underestimates the self-diffusion coefficient value. However, the minima of the curves for each species of the mixtures are less pronounced than in the experiment. Moreover, the crossing point where the self-diffusion coefficient of water and methanol molecules are equal is predicted by the simulation at a high X m (X m ≈ 0.6), whereas in the experiment this crossing occurs at an even higher methanol fraction, around 0.7. It is not straightforward to establish any correlation between the behavior of mixing properties and the self-diffusion coefficients of the species. Experimentally, the minimum values of excess mixing volume and enthalpy, as well as of self-diffusion coefficients, occur at a slightly prevailing water fraction, cf. figures 2, 3 and 6. As concerns the excess volume, visually it looks as if the minimum position was at X m = 0.5. On the other hand, simulations involving the OPLS/UA methanol model indicate that the minima for ∆H mix and D w are located at a prevailing methanol fraction. In order to explore this discrepancy more in detail we resort to the notion of hydrogen bonding.
The structural, mixing and dynamic properties in hydrogen bonded liquids are all certainly locally connected with the characteristics of the H-bonds, and globally with the nature of the hydrogen bonded network that is formed in a given system. As a step towards understanding this connection, we have calculated characteristics related to H-bonding. We decided to use a rather popular geometric criterion applicable to water and methanol that involves two distances and one angle, see e.g., [37,38] . Finally, the angle O x − H x · · · O y should be smaller than the threshold value, see e.g., references [39][40][41] concerning the analysis of the structure of liquid methanol. If the coordinates of two molecules fulfill the above conditions, they are considered to be H-bonded. The values used in this work are the natural ones coming from the first (intermolecular) minima of the corresponding radial distribution functions for atoms O-O and O-H belonging to different species at each calculated composition. There is a certain freedom concerning the choice of a threshold value for the angle. We have chosen to use the most popular threshold value for the angle between two vectors, r O x O y r O x H x , equal to 30 • that yields, as we will see below, a higher fraction of H-bonded molecules in comparison with the equally successful, more restrictive criterion for the angle between vectors, r H x O x r H x O y , larger than 150 • , see the angle θ in figure 1 of reference [37].
In the two panels of figure 7, the average number of hydrogen bonds per molecule is given as a function of composition for two different geometric definitions of H-bonds. Both definitions lead to similar trends whereas the actual numbers at each composition are a bit different. The values are normalized per the number of molecules of the first species marked next to a given line (wwater; mmethanol).
That is, for instance, the line "w-m" gives the number of hydrogen bonds, 〈n〉 HB , for a water molecule that forms connections with methanol ones (normalized by the number of water molecules). In figure 7 (b), the values for the UA and AA models of methanol are also compared: correspondence between the two types of potentials is nearly perfect in this respect. In both panels of figure 7 the vertical dashed line marks an approximate composition (X m ≈ 0.60) at which the excess mixing volume is at minimum, cf. "w-w" decays in the entire range of composition reflecting the shrinking of the system due to breaking of H-bonds, whereas the "w-m" line shows an increasing average number of bonds and leading to "expansion" of the system due to an increasing number of bonds. The "m-m" and "m-w" lines exhibit similar trends. Consequently, in two regions where one of the two average numbers of bonds prevail, the excess mixing volume decreases or increases. In the intermediate interval of compositions, where the average numbers of bonds are balanced, ∆V mix exhibits its minimum. These trends are independent of the definition of the H-bond. If ∆V mix is in the interval of its minimum value, the self-diffusion coefficients are also approximately around their respective minima. This picture comes from computer simulations of both of the two water models, combined with OPLS/UA methanol; we show only the case of SPC/E to avoid an unnecessary repetition. It is worth mentioning that the excess mixing enthalpy is not straightforward to incorporate into the interpretation of the above data. The reason is that ∆H mix has "geometric" and energetic terms. Therefore, combining them to interpret the observed trends using a purely geometric definition of H-bonds does not seem to be appropriate.
One interesting conjecture coming from the results in figure 7 is that as a general tendency, both water and methanol molecules prefer to coordinate water molecules via H-bonding, similarly to what was found for the mixtures with AA methanol molecules [1]. This can be discerned from the non-linear dependence of the number of H-bonds per molecule on the methanol concentration. In the non-preferential case, the number of both water and methanol H-bonded neighbors would change linearly with concentration. The deviation from linearity is positive for the "w-w" and "m-w" curves, i.e., there are more water molecules hydrogen bonded to both water and methanol molecules than it would be proportional to the composition. By contrast, the deviation is negative for the "w-m" and "m-m" curves, i.e., the number of H-bonded methanol molecules to both water and methanol molecules is smaller than it would be proportional to the methanol concentration. Interestingly, the effect is more pronounced for the less strict definition of H-bonds, as exemplified by figure 7 (a). By a close inspection of the fractions of methanol and water molecules in different bonding states, as a function of composition (see figure 8), a very similar conclusion may be drawn. These graphs provide the ratio of water or methanol molecules that have a specific number of H-bonded water or methanol neighbors at a given composition. Water can form a maximum of four hydrogen bonds, whereas methanol can form three. If we compare water-water ["w-w" in figure 8 (a)] and water-methanol ["w-m" in figure 8 (c)] curves, the following observations can be made. As concerns panel (a) of figure 8 that describes the bonding solely within the water subsystem of the mixture, the following trends are worth mentioning. Namely, when the methanol concentration increases (starting from pure water up to pure methanol), the fractions of water molecules with four and three "w-w" bonds monotonously decrease whereas the fraction of molecules that do not participate in bonding substantially grows. Only the fractions of doubly and singly bonded water molecules exhibit maxima at a respective composition, X m . In other words, clusters (and branched structures) of exclusively water molecules transform into chains, and thus dimers as well as "free" molecules become abundant. This is a crude description of how a water subsystem of a mixture shrinks. Actually, the minimum value of the excess mixing volume is located in between the two maxima observed for doubly and singly bonded waters.
On the other hand, see figure 8 (c), the fractions describing the bonding of water molecules solely

43602-10
with methanol species exhibit the following trends. The average number of waters initially free of "w-m" bonds rapidly decreases at the expense of the formation of bonded structures. Only the fractions of water molecules with a single bond and double bonds with methanols pass through a maximum value, other fractions with higher number of bonds grow monotonously. Thus, the expansion of the water subsystem occurs via substitution of water molecules by methanols in a shell of neighbors of a given water molecule.
In close similarity to what we have seen in panel (a), the minimum of the excess mixing volume occurs at X m located in between two maxima of the curves describing singly and doubly bonded molecules. Similar interpretation can be developed to describe the changes of the bonding within the methanol subsystem and for the changes of the fraction of "m-w" bonds, panels (b) and (d) of figure 8.
Our final focus in the present study is in the relaxation properties of molecules in mixtures. Those certainly are also influenced by H-bonding discussed above. In close similarity to the previous studies [15,42], we would like to evaluate the reorientational correlation functions determined from where P l is the lth order Legendre polynomial and e α is the unit vector which points along the α axis in the molecular reference frame. The average is performed over all molecules belonging to a particular species. A typical plot for a mixture of a given composition describing the dependence of the autocorrelation function on time is given in figure 9. In the present calculations we used the axis parallel to the molecular dipole moment µ which is related to dielectric relaxation. Moreover, we calculated C µ l for l = 2 for the two species of the mixture. In addition, we performed calculations by using the O-H vector (H is the hydrogen atom belonging to the hydroxyl group of the molecule) as α axis to obtain the corresponding relaxation time. The single-molecule reorientational times follow from integration of the autocorrelation function defined above: To perform the integration we split the time interval into two pieces, actually it is a commonly accepted procedure [10,42]. Over a short-time interval, the integration is performed numerically; the upper limit of this interval is chosen between 4 and 5 ns approximately. For a long-time interval, the logarithm of the autocorrelation function is plotted: it behaves as a straight line permitting to analytically evaluate the relevant contribution [42].
Single molecule reorientational times, as calculated for both kinds of the axes mentioned above, are reported in figure 10, as a function of composition. Two trends are qualitatively similar in the behavior of  It is difficult to establish how accurate all these results are due to a lack of well established experimental data. We included some points obtained by Ludwig [43] in figure 10 (b). At present it seems impossible to perform a comprehensive validation of different models versus such kind of experimental data.

Summary
To summarize this study, we have used molecular dynamics simulations in the NPT ensemble to study some thermodynamic properties, as well as the dynamic and dielectric properties of water-methanol mixtures at room temperature and at ambient pressure in the entire range of their composition. The SPC/E and TIP4P/Ew models combined with the OPLS/UA potential model for methanol are used. This is done as a first step of systematic studies of mixtures of water and a set of alcohols using nonpolarizable models. Our principal findings concern the behavior of density, excess mixing volume and excess mixing enthalpy on composition. Also, we have evaluated the self-diffusion coefficients of both species, reorientational relaxation times and trends of behavior of hydrogen bonding between molecules. Whenever it seemed possible, we looked for a unified interpretation of different properties and performed comparisons with available experimental results. It would be interesting to establish a relation between the structural properties and hydrogen bonding on the one hand and the dielectric constant on the other hand, possibly in the spirit of previous studies of diffusion and power spectra, see e.g., [44]. We would like to note that in our opinion an ample set of properties of the model mixtures in question is required to establish if the modelling is consistent (with experimental data) and successful. A combination of experimental data, computer simulation results and reverse Monte Carlo modelling along the lines developed in [45,46] may be profitable in this respect. Possible modifications of parameters and combination rules nevertheless should be explored more in detail. At present, it seems that each specific parametrization leads to a successful description of a particular property, such state of art being unsatisfactory. A more sophisticated modelling of alcohols, e.g., OPLS/AA, and the use of possibly more refined water models is necessary to extend the present study and at least make stronger conclusions about various thermodynamic properties of these mixtures. Unfortunately, it seems that some of them are intrinsically impossible to describe very well within the framework of non-polarizable models.