Article pubs.acs.org/JPCB

Investigation of the Structure of Ethanol−Water Mixtures by Molecular Dynamics Simulation I: Analyses Concerning the Hydrogen-Bonded Pairs Orsolya Gereben* and László Pusztai Research Institute for Solid State Physics and Optics, Hungarian Academy of Sciences, P.O. Box 49, H-1525 Budapest, Hungary S Supporting Information *

ABSTRACT: Series of molecular dynamics simulations for ethanol−water mixtures with 20−80 mol % ethanol content, pure ethanol, and water were performed. In each mixture, for ethanol the OPLS force field was used, combined with three different water force fields, the SPC/E, the TIP4P-2005, and the SWM4-DP. Water potential models were distinguished on the basis of deviations between calculated and measured total scattering X-ray structure factors aided by ethanol−water pair binding energy comparison. No single water force field could provide the best agreement with experimental data at all concentrations: at the ethanol content of 80% the SWM-DP, for 60 mol % the SWM4-DP and the TIP4P-2005, whereas for the 40 and 20 mol % mixtures TIP4P-2005 water force field provided the closest match. Coordination numbers and hydrogen bonds/molecule values were calculated, revealing that the oxygen−oxygen first coordination numbers strongly overestimate the average number of hydrogen bonds/molecule. The center-of-molecule distributions indicate that the ethanol−ethanol first coordination sphere expands with increasing water concentration while the size of the first water−water coordination sphere does not change. Various two and three-dimensional distributions were calculated that reveal the differences between simulations with different water force fields. Detailed conformational analyses of the hydrogen−bonded pairs were performed; drawings of the characteristic molecular arrangements are provided. tried so far. The simpler so-called “united atom” approach is when not all the atoms are represented; (usually) CH3 and CH2 groups are combined in one site, leaving only four interaction sites in total for an ethanol molecule. For instance, the AMBER united atom force field13 falls in this category, as well as the original version of the frequently used Optimized Potentials for Liquid Simulations (OPLS) parameters, developed by Jorgensen et al.,14 that have been applied, for example, in the work of Wensink et al.12 OPLS formed the basis of further optimized potential parameters for specific systems.15,16 However, the united site models, although computationally less time-consuming, have disadvantages; as not all the atoms are represented, they miss some of the partial radial distribution functions and stereochemical information, like some bond angles and dihedral angles. Therefore, with increasing computer capacity, development of all-atom force fields started; a well-known specimen is the OPLS-AA (All Atom) force field for several types of molecules.17,18 A further step toward increasing the accuracy of the models is the inclusion of polarization effects. The Polarizable Intermolecular Potential Function (PIPF) force field has been developed by Gao et al.19 and significant polarization effects were found for pure ethanol. Noskov et al. proposed another polarizable

1. INTRODUCTION Aliphatic alcohols have been targets of scientific investigations for a long time; for instance, they were among the first to be studied by X-ray diffraction.1−3 They have an amphiphilic character, consisting of both hydrophobic and hydrophilic parts, which duality makes their behavior rather complex. Their molecules can form hydrogen bonds and therefore they may play important roles in living organisms as well, and this is also a reason for the extensive investigations mentioned above. Focusing on the structure, analyses now include various computer simulation techniques beside diffraction experiments, for example, molecular dynamics (MD) simulations,4−6 configurational-bias Monte Carlo simulation,7 and reverse Monte Carlo (RMC) modeling.8 These approaches have been aimed at understanding the atomic level structure, including characteristics beyond pair correlations, and aggregation7 of these materials. Ethanol is used widely in different branches of industry and in everyday life, as well, as a chemical reagent, solvent, paint stripper, fuel, and a component in alcoholic beverages. Pure ethanol and its mixtures with water have therefore been investigated by several experimental techniques, for example, NMR spectroscopy9 and dielectric relaxation.10 Molecular dynamics (MD) computer simulations are very popular among the simulation methods studying ethanol and its water mixtures,11,12 even though choosing appropriate potential functions is not easy. Several approaches have been © XXXX American Chemical Society

Received: October 18, 2014 Revised: January 15, 2015

A

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B ethanol model,20 based on a Drude oscillator, which provided dielectric constant for mixtures in good agreement with experimental data. For water molecules in alcohol−water mixtures, various rigid water potentials have been tested. The SPC/E force field of Berendsen et al.21 was used in the work of Tarek et al. for the simulation of a 0.1 M ethanol−water solution,22 the TIP4P force field of Jorgensen et al.23 was utilized in the work of Zhang4 and Wensink.12 The TIP3P23 force field was made use of in the work of Zhang et al.6 Noskov et al.20 applied the SWM4-DP force field in their work. In a wider context, ethanol−water mixtures have been investigated from different viewpoints; apart from structural studies, dynamic properties were also calculated. The properties of the air/solution interface were examined by Tarek et al.22 for a 0.1 M ethanol−water mixture; they found that the number of water molecules involved in hydrogen bonding with ethanol molecules on the surface decreased by a factor of 2 compared to the bulk. Wensink et al.12 in their study of ethanol−water mixtures concentrated on the calculation of densities, energies, shear viscosity, and diffusion properties; they found that although for all the properties examined, the excess values are reproduced qualitatively, although all of them are underestimated. They suggested that although the applied OPLS and TIP4P potentials describe well the pure substances, they fail to reproduce properties of the mixtures quantitatively. The complex dielectric permittivity at eight different temperatures and five different concentrations was investigated by Petong et al.;10 they found that the assumption of two relaxation regions can describe the behavior of these mixtures. Noskov et al.20 studied the hydrophobic hydration using a polarizable ethanol model, combined with the SWM4-DP24 water model, and found that both ethanol and water prefers water in their first coordination shell. In the present study, we performed a series of molecular dynamics simulations for ethanol−water mixtures of 20, 40, 60, and 80 mol % ethanol content and for pure ethanol and water as reference. The OPLS-AA force field18 was used to describe the ethanol molecule. It was combined with three different water force fields: (1) the simple three-site SPC/E21 force field; (2) from the family of four interaction site (Bernal-Fowler25 geometry) models the TIP4P-200526 one was selected, as it is an improvement considering several physical properties26 over the original TIP4P potential;23 (3) in order to take polarizability into account, the SWM4-DP force field24 was chosen. All the three water potentials were used for all mixtures of this study. Our first aim was to find the most appropriate force field combination in the entire concentration range, or, should this not be possible, for each concentration separately, mainly based on a direct comparison of the calculated (from MD simulation) and measured total scattering X-ray structure factor (TSSF), F(Q). To a lesser extent, a comparison of our MD based ethanol−water dimer binding energies with related results of ab initio quantum chemical calculations27 was also taken into account. This paper reports the comparison of the TSSFs, along with some of the partial radial distribution functions (PRDF), g(r), in Section 3.2, as well as the comparison of the ethanol−water dimer binding energies (in Section 3.3). Coordination numbers and hydrogen bond statistics are provided in Sections 3.4 and 3.5, respectively. center-of-molecule distributions are discussed in Section 3.6, and finally in Section 3.7 a conformational analysis of hydrogen-bonded pairs is given for the most

successful models. A detailed study of the hydrogen-bonded network, identification of cyclic and acyclic parts, and complete ring analyses will be reported separately.

2. METHODS 2.1. Molecular Dynamics Simulations. For the MD simulation, the GROMACS 4.0 simulation package28 has been used throughout for NVT ensembles at T = 293 K with time step dt = 2 fs. The equations of motions were integrated via the leapfrog algorithm. The particle-mesh Ewald29,30 (PME) algorithm was applied in order to account for long-range electrostatic interactions. The cutoff value of 11 Å used by Jorgensen et al. in the parametrization process for alcohols18 were adopted. The temperature was maintained by a Nose− Hoover thermostat31 with τT = 0.5 ps. The geometry of water molecules was maintained by the SETTLE algorithm.32 Flexible ethanol molecules were used, including all the possible dihedral angles. The simulation length was 2000 ps in each case. Particle configurations in the “production” phase were collected 20 ps apart and in the end, 76 configurations were used for calculating the average radial distribution functions and structure factor. Cubic simulation boxes have been applied; the number of molecules (determined form the appropriate solution densities) is given in Table 1, along with the number densities and solution densities for all the studied concentrations. Table 1. Number of Ethanol (nE) and Water (nW) Molecules and Number Densities for the Different Ethanol−Water Mixtures.a Et100 Et80 Et60 Et40 Et20 Et0

nE

nW

ρE (Å−3)

ρW (Å−3)

ρ (Å−3)

ρliq (g/cm3)

1331 1226 1093 889 571 0

0 315 714 1326 2280 3993

0.0103 0.0096 0.0088 0.0073 0.0047 0.0000

0.0000 0.0025 0.0057 0.0109 0.0188 0.0330

0.0930 0.0940 0.0960 0.0985 0.0990 0.0990

0.790 0.799 0.816 0.837 0.841 0.987

ρE and ρW are the molecular number densities of ethanol and water molecules, respectively, whereas ρ denotes the total atomic number density of all the atoms in the simulation cell. The density of the liquids also given in g/cm3. a

For the MD calculations, all-atom force field parameters were needed in order to make the calculation of all the PRDF-s and the proper total scattering X-ray [F(Q)] structure factor possible. For ethanol, the OPLS-AA18 force field was used; concerning water, three different water force fields, SPC/E,21 TIP4P-2005,26 and SWM4-DP,24 were applied at each concentration. All in all, 12 ethanol−water mixtures, 3 pure water, and 1 pure ethanol simulations were carried out. (The pure ethanol and water simulations were needed as references.) The equilibrium molecular geometries for the three water force fields are provided in Table 2, whereas full specification of the OPLS-AA force field parameters can be found in the literature.18,28 The nonbonded interaction parameters between ethanol and water sites were calculated as σ12 = √(σ1σ2), and ε12 = √(ε1ε2) (combination rule “3”, according to the GROMACS convention28). Calculations reported here will be identified by their ethanol content and the force field applied, so, for example “Et80 SPC/ E” will refer to the 80 mol % ethanol containing mixture where the MD simulation was performed using the SPC/E force field. B

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

3. RESULTS AND DISCUSSION 3.1. Water. Although detailed investigation of pure water is not the main subject here, the differences observed between the different ethanol−water mixtures in terms of the watercontaining partials made performing MD simulations necessary with our chosen three force fields for pure water, as well. In this subsection, some particulars of the total and partial structure factors and PRDF-s are given briefly. It may be in order to note in advance that concerning the structure of pure water at ambient conditions, the SWM4-DP and TIP4P-2005 potentials have proved to be significantly more consistent with neutron39 and X-ray40 diffraction experiments than the SPC/E one. Concerning r-space information, there are only negligible differences between the O−O partial RDFs of the three MD simulations, that manifest in the increasing height (in order SPC/E → TIP4P-2005→ SWM4-DP) of the first (pre)peak of the structure factor around 2.1 Å−1 (not shown). The intramolecular parameters of the SPC/E force field are different from the SWM4-DP and TIP4P-2005 force fields (see Table 2), which the latter two are identical. This is reflected by the first peak positions of the O−H and H−H PRDFs. There is a slight difference in terms of the intermolecular O−H structure between the SPC/E and the other two models, reflected by the difference of the second peak of the PRDF (not shown). The larger differences between the SPC/E and the other two water potentials in terms of the structure factors (see Figure 1) are coming mainly from the different intramolecular structure; no visible differences between the F(Q) of SWM4-DP and TIP4P-2005 can be observed. While there are only minor differences in terms of the intermolecular H−H PRDF (not shown), there are visible discrepancies even between the corresponding SWM4-DP and TIP4P-2005 force fieldgenerated partial structure factors (see Figure 1a). The total structure factors are shown in Figure 1b (standard deviations are between 0.005 and 0.06, not shown in figure, as they would impair visibility); the main difference can be seen around the first peak. In summary, it may be concluded that the main difference of the TSSF-s originate mainly from the intramolecular and to a smaller extent from the intermolecular H− H interactions. 3.2. Comparison of X-ray Structure Factors and Partial Radial Distribution Functions of Water−Ethanol Mixtures. Experimental and MD-simulated structure factors for mixtures of 80, 60, 40, and 20% ethanol content are displayed in Figure 2. (Standard deviations for the simulated F(Q) data points for the ethanol-containing samples varied between 0.01 and 0.1; the error bars are not shown in the figures, as they would impair visibility.) Before going into details, it may be stated that the overall “visual” consistency

Table 2. Intramolecular Equilibrium Distance Parameters for the Various Water Models Used in the Present MD Simulations d0 (Å) O−H H−H

SPC/E

SWM4-DP

TIP4P-2005

1.0000 1.6330

0.9572 1.5139

0.9572 1.5139

2.2. Structure Analysis and Binding Energies. The PRDFs and the TSSF were calculated by the g_rdf_rmc software, which is our modified version of the g_rdf program provided in the GROMACS package.28 Software used for the calculation of dimer binding energies, various distributions and statistics concerning the hydrogen bonding (see below) has been developed specifically for the present work. Note that an “absolute” definition of hydrogen bond is not available; here, H-bonds have been defined similarly as in our previous work.33 The LOOSE definition is based solely on geometric considerations: the H-bonding O---O and O---H distances had to be in distance ranges specified by the appropriate first intermolecular coordination spheres (i.e., the neighborhood of the first intermolecular maxima). The upper limiting values were defined by the first minima of the O−O, and by the second (i.e., first intermolecular) minima of the O− H PRDFs. The maximum and minimum positions have changed somewhat from simulation to simulation: on average, the “H-bond” ranges were set to be between 2.4 and 3.6 Å for O−O and 1.4−2.7 Å for O−H pairs (see the Supporting Information Tables 1 and 2 for details). As the origin of hydrogen bonding is considered mainly electrostatic, the strength (and stability) of the H-bond does not only depend on the O---H distance, but on the O−H---O angle, as well. Therefore, a STRICT condition, similarly to Chen and Siepmann,7 has also been applied, where only O−H---O angles larger than 120° were considered to form a hydrogen bond. The 120° threshold was determined based on actual O−H---O angle distributions; it also corresponds to values previously used in earlier works.7,33 2.3. Experimental Data from X-ray Diffraction. For a thorough comparison with experiment, the X-ray diffraction data measured by Temleitner34 in the SPring-8 Synchrotron Facility (Hyogo, Japan), using the BL04B2 high-energy liquid and amorphous X-ray diffraction beamline,35 were taken for the ethanol−water mixtures. Here, experimental data in the form of total scattering structure factors were made use of after the raw data had been corrected for background and Compton scattering via standard procedures.36,37 X-ray diffraction data for pure ethanol8 and pure water38 were published earlier.

Figure 1. (a) H−H partial structure factors and (b) the total experimental and calculated X-ray structure factor of the SPC/E, TIP4P-2005, and SWM4-DP water simulations. C

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 2. Experimental and MD simulated total structure factors of the SPC/E, TIP4P-2005, and SWM4-DP MD simulations for the (a) Et80, (b) Et60, (c) Et40, and (d) Et20 ethanol−water mixtures.

signal of an unrealistically high intensity, as well as larger deviations from experiment than those corresponding to the other two water potentials. Comparing the water−water oxygen g(r), OW−OW, as a function of concentration (Figure 3a) for the same force field, the position of the first maximum does not change, although its magnitude increases. This is analogous to the MD simulation data of Zhang et al.,4 and to the experiment-based findings of Dixit et al.41 who considered a 70−30% methanol−water mixture. The position of the second maximum shifts gradually from ∼4.5 (pure water) to ∼4.9 Å (20% water) with the increasing ethanol concentration, while its height also increases somewhat, contrary to the findings of Dixit41 for their methanol−water mixture, where both the height and the position remained basically unchanged. This difference between the two alcohols might be attributed to the larger size of the ethanol molecule, involving a larger hydrophobic part, thus disrupting the Hbonding network of water to a larger extent. The difference is probably not the effect of a decreased ability to form H-bonds, as both methanol and ethanol can form the same number of Hbonds per molecule. Comparing the behavior of the ethanol−water OE−OW PRDFs, the same tendency can be observed for the simulations with the same water force field and different concentrations; the position of the second maximum is shifting toward larger distances with increasing ethanol concentration, although in a less regular fashion as compared to the OW−OW case (see Figure 3b). In the figure, the OE−OE PRDF of pure ethanol and the OW−OW PRDF of pure TIP4P-2005 water are also given for comparison. The position of the first OE−OE maximum is at 2.79 Å according to our calculation, which is in better agreement with the value of 2.8 Å of Narten et al.42 based on X-ray diffraction data than with that of 2.74 Å found by Jorgensen14 by molecular dynamics simulation. The second O−O maximum in pure water is around 4.45 Å while it can be found at 5.05 Å in pure ethanol. Interestingly, the second OE− OW peaks of the Et20, Et40, and Et60 mixtures hardly differ

between measurement and simulation is (perhaps surprisingly) of high level. Discrepancies between simulations performed with various water force fields at the same concentration, as well as between simulation and experimental data are increasing with the increasing water concentration. Interestingly, differences between simulations for pure water are smaller than differences between Et20 simulations; the three water force fields applied here interact differently with the ethanol force field and water− ethanol interactions appear to magnify differences between water models. The success of the different water force field simulations may be characterized by calculating the squared difference, χ2, between the MD (averaged over many time frames) and experimental structure factors; the values are given in Table 3. Table 3. Goodness-of-Fit, χ2, for the Different MD Simulations of Ethanol−Water Mixtures and Pure Ethanol and Watera Et100 Et80 Et60 Et40 Et20 Et0

χ2 SPC/E

χ2 TIP4P-2005

χ2 SWM4-DP

2.0 1.6 1.5 1.5 1.6 1.7

2.0 1.5 1.3 1.1 1.3 1.3

2.0 1.2 1.3 3.2 3.9 4.7

a The MD simulation with the smallest χ2 for each concentration is highlighted with bold.

(Note that it would be hard to make such a distinction on the basis of mere visual inspection.) For the Et80 calculation, the SWM4-DP force field produced the smallest χ2, whereas for the Et60 simulation both the SWM4-DP and TIP4P-2005 force fields gave equally good results. At higher water content, the best agreement was achieved by using the TIP4P-2005 force field. It seems that at lower water ratios the polarizability of the SWM4-DP force field has beneficial effects; on the other hand, at higher water concentrations it causes a small angle scattering D

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. (a) OW−OW and (b) OE−OW g(r) partials for the TIP4P-2005 simulations in comparison with the OW−OW partial of pure TIP4P-2005 water and with the OE−OE partial of pure ethanol (only in case of b)). Series names are ordered in the order of decreasing first peak height.

Table 4. (a) The Binding Energies (with Standard Deviations, in kJ/mol) of the Hydrogen-Bonded Ethanol−Water Pairs for the MD Simulations Performed in This Work, for All the Mixtures (see the description of the abbreviations in the text) and (b) The Binding Energies (in kJ/mol) of the Ethanol−Water Dimers from the Literature a)

SPC/E

LOOSE

TIP4P-2005

SWM4-DP

SPC/E

TIP4P-2005

EHA

SWM4-DP

EHD

Et80 Et60 Et40 Et20 LOOSE

22.34 21.87 21.44 21.34

± ± ± ±

0.36 0.26 0.28 0.33

24.54 ± 0.39 23.96 ± 0.29 24.16 ± 0.22 23.78 ± 0.26 EgHA

24.20 25.63 27.15 27.19

± ± ± ±

0.31 0.28 0.10 0.13

19.09 18.67 18.11 17.98

± ± ± ±

0.47 0.40 0.45 0.45

17.10 ± 0.44 16.70 ± 0.34 17.22 ± 0.30 17.09 ± 0.31 EgHD

18.92 20.19 22.24 21.91

± ± ± ±

0.08 0.27 0.22 0.26

Et80 Et60 Et40 Et20 STRICT

21.90 21.33 20.77 20.69

± ± ± ±

0.38 0.31 0.42 0.31

23.85 ± 0.37 23.37 ± 0.33 23.53 ± 0.28 23.06 ± 0.26 EHA

23.60 24.83 26.42 26.04

± ± ± ±

0.40 0.34 0.16 0.19

18.93 18.40 17.89 17.64

± ± ± ±

0.48 0.42 0.42 0.41

17.06 ± 0.41 16.79 ± 0.31 17.22 ± 0.31 16.92 ± 0.40 EHD

18.72 19.77 21.08 21.30

± ± ± ±

0.49 0.35 0.24 0.28

Et80 Et60 Et40 Et20 STRICT

22.50 22.07 21.68 21.47

± ± ± ±

0.36 0.24 0.27 0.33

24.62 ± 0.37 24.07 ± 0.29 24.17 ± 0.22 23.80 ± 0.26 EgHA

24.30 25.72 27.15 27.19

± ± ± ±

0.28 0.26 0.10 0.13

19.20 18.82 18.31 18.09

± ± ± ±

0.48 0.39 0.40 0.46

17.15 ± 0.41 16.80 ± 0.35 17.23 ± 0.29 17.10 ± 0.30 EgHD

19.08 20.25 22.27 21.94

± ± ± ±

0.44 0.25 0.22 0.24

Et80 Et60 Et40 Et20

22.02 21.52 21.01 20.81

± ± ± ±

0.37 0.28 0.31 0.30

23.95 23.51 23.54 23.09

± 0.34 ± 0.32 ± 0.27 ± 0.26 EHA

23.68 24.89 26.42 26.05

± ± ± ±

0.40 0.34 0.16 0.19 EHD

19.00 18.54 18.09 17.74

± ± ± ±

0.48 0.42 0.38 0.40

17.11 16.87 17.23 16.94

0.48 0.34 0.23 0.27

EgHA

18.86 ± 19.83 ± 21.83 ± 21.36 ± EgHD

23.64 23.51 20.21 19.58

18.83 18.83 14.48 14.69

b) MP427 CCSD(T)27 CPMD-BLYP44 ADF-BLYP44

22.51 23.38 19.58 18.87

17.87 17.91 16.44 15.86

± ± ± ±

0.40 0.30 0.31 0.39

3.3. Hydrogen Bonded Ethanol−Water Pair Binding Energies. The binding energies of the hydrogen bonded ethanol−water pairs were also calculated for all the ethanol− water mixture simulations. As our ongoing cluster analyses of these systems (to be reported separately) revealed, there are hardly any dimers (clusters of size two) in the mixtures, so these hydrogen-bonded pairs are mostly parts of larger clusters and not true dimers; for the sake of simplicity, we will refer to them as dimers from now on. It has to be noted that ethanol can exist in two stable conformations, trans or gauche. On the basis of the microwave study of Kakar et al.43 and on theoretical calculations,27,44 the trans conformer is slightly more stable. Taking also into account whether water or ethanol is the hydrogen donor, four different kinds of dimers can be distinguished: the trans ethanol−water, where ethanol is either H acceptor (EHA), or H donor (EHD), and similarly, the gauche ethanol−water,

from each other, although they are at a larger distance than the OW−OW second peak of pure TIP4P-2005 water. (The same tendency can be observed for the OE−OE g(r) partials of the mixtures, not shown.) Only the Et80 mixture, containing the smallest amount of water, shows a shift to an even larger distance, approaching the maximum position found in pure ethanol. This trend can be explained as follows: as long as there is a sufficient amount of water in the system (>20%), the Hbonding structure of water dominates, keeping the ethanol and water oxygen atoms, and even the ethanol−ethanol oxygen atoms, closer to each other. When the ethanol concentration reaches a threshold, the ethanol H-bonding network starts to dominate; the molecules move away from each other, closer to the equilibrium value of the pure ethanol. In order to fully understand these phenomena, a very detailed investigation of the H-bonding networks of these mixtures is presently being carried out. E

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

are more stable than the ones where ethanol is the donor, similarly to the works of both Fileti27 and van Erp.44 In our case, the trans ethanol conformer has larger binding energy than the gauche one, regardless of whether ethanol is H donor or acceptor. The squared difference, χ2, of the binding energies for each force field and concentration value, in comparison with the appropriate MP4 and CCSD(T) results of Fileti,27 were calculated. It depends on the type of the dimer and on the concentration which force field produces the value closest to the ab initio results. Calculating the partial sum for an ethanol concentration over the four dimer types for each force field (see Table 5) produced similar results as the analysis of the X-ray

where ethanol is H acceptor (EgHA) or H donor (EgHD). Ethanol molecules were considered to be trans if the CCOH dihedral angle was between 120 and 240°; otherwise, the molecule was categorized as gauche. The average dimer binding energies in kJ/mol are given in Table 4a for all the four concentrations and three water force fields and for both the LOOSE and STRICT definitions of the H-bond. The following observations can be made from studying the binding energy values: (1) The binding energy shows no significant dependence on the concentration. For the SPC/E force field there is a slight increase of the binding energy for all the four kinds of dimers with increasing ethanol concentration, while for SWM4-DP more or less a reverse tendency can be observed. For TIP4P2005, there is no tendency at all. (2) Comparing the values of the LOOSE and STRICT HB conditions for the same concentration and force field, the binding energy is almost always larger for the STRICT condition. This can be explained by the fact that in case of the STRICT condition, more dimers have a geometry closer to the ADF-BLYP ab initio44 optimized H-bond angle (varying between 168 and 174°), depending on the type of the dimer. Around 29% of the pairs for Et80 and 24% for Et20 have close to optimal geometry in case of the LOOSE condition. No experimental values for the dimer binding energies could be found in the literature. The energies of the ethanol−water dimers were investigated by Masella et al.45 applying MP2 level ab initio calculations, by van Erp et al.44 by DFT-based Car− Parrinello MD and ADF-BLYP simulations, and by Fileti et al.27 using a range of different level ab initio calculations. Unfortunately, the results obtained by different methods and authors are not much in agreement, but some trends can be revealed. On the basis of the work of Fileti27 and van Erp,44 the dimers where the ethanol is the H-acceptor are more stable than the ones where the ethanol is the H-donor, regardless the conformation of the ethanol molecule (see Table 4a). On the other hand, the earlier work of Masella45 found the ethanol hydrogen donor complex more stable. This discrepancy was attributed to the limited basis set in the MP2 calculation used by Masella45, and therefore their results are not considered here any further. Considering the trans and gauche conformer, in both Fileti’s and van Erp’s work EgHA was more stable than EHA. In the work of Fileti,27 EgHD was more stable than EHD, while the reverse order was found by van Erp.44 The binding energies (see Table 4b) coming from the CPMD-BLYP and ADF-BLYP calculations of van Erp44 have significantly lower values than the values calculated by Fileti.27 By the authors of that work, this was attributed to the inability of the gradient-corrected BLYP functionals to account for the dispersion forces. Therefore, the binding energies calculated in the present work will be compared to the highest level ab initio results of Fileti.27 However, it has to be kept in mind that while the ab initio results are calculated for optimized geometry gas phase dimers, our results are average liquid phase values of ethanol− water mixtures with varying concentration and calculated for hydrogen-bonded pairs of mainly larger clusters. It is not even certain whether the equilibrium geometry of the dimers in the liquid phase should be the same as in the gas phase. Bearing all this in mind, the binding energies calculated in this work are in very good agreement with the MP4 and CCSD(T) results of Fileti.27 For all the three water force fields and all concentrations, the pairs where ethanol is the acceptor

Table 5. Subtotal Squared Difference, χ2, of the Binding Energies of the H-Bonded Ethanol−Water Pairs Compared To Ab Initio Values of Fileti.27,a χ2 ab initio level − HB condition

ethanol conc.

CCSD(T) − MD LOOSE HB

Et80 Et60 Et40 Et20 sum Et80 Et60 Et40 Et20 sum Et80 Et60 Et40 Et20 sum Et80 Et60 Et40 Et20 sum

CCSD(T) − MD STRICT HB

MP4 − MD LOOSE HB

MP4 − MD STRICT HB

SPC/E

TIP4P2005

SWM4-DP

5.1 7.8 12.2 13.5 38.6 4.7 6.6 9.8 12.2 33.3 4.6 6.6 10.3 11.5 32.9 4.4 5.7 8.3 10.3 28.8

5.3 6.0 3.7 4.7 19.6 5.3 5.5 3.6 4.6 19.0 7.9 7.7 5.7 6.2 27.6 8.0 7.4 5.7 6.1 27.3

1.7 12.9 46.5 43.0 104.1 2.2 13.9 50.7 43.6 112.6 4.0 17.4 53.4 50.1 124.9 4.7 18.5 57.6 50.7 133.8

a

Summation is over the dimer types for a given ethanol concentration, and their total over the concentration range for the three applied water force fields is also shown. The smallest values are highlighted by bold characters.

structure factors. Namely, the Et80 SWM4-DP case gave binding energies closest to the CCSD(T) (both for the LOOSE and STRICT HB conditions), and to the MP4 values when using the LOOSE HB condition. In the case of the MP4−Et80 MD STRICT HB condition, SPC/E gave the closest binding energy. In the case of CCSD(T), both for the LOOSE and STRICT HB conditions, TIP4P-2005 proved to be the best for Et60-Et20, similarly to the structure factor results. Comparing to MP4, SPC/E proved to be the best for Et60, and TIP4P2005 for Et40-Et20, for both HB conditions. Summing the χ2 values for the same dimer over the concentration range for each force field revealed that for the CCSD(T) ab initio data, TIP4P-2005 produced the smallest χ2 for the EHA and EgHA, and SPC/E for the EHD and EgHDtype dimers (for both the LOOSE and STRICT HB conditions). Comparing to the MP4 values, for the EHA pairs the χ2 of the SPC/E force field was the smallest for both F

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 4. Upper panels: Oxygen−oxygen coordination numbers for the first coordination spheres for (a) ethanol−ethanol, ethanol−water and ethanol−(ethanol + water), denoted as OE−OE, OE−OW, and OE−OX plotted against the ethanol number density and for (b) water−water, water− ethanol, and water−(ethanol + water), denoted as OW−OW, OW−OE, and OW−OX, plotted against the water number density. Lower panels: similarly the oxygen−hydrogen coordination numbers for the first coordination spheres for (c) ethanol−ethanol, ethanol−water, ethanol−(ethanol + water), denoted as OE−HE, OE−HW, and OE−HX plotted against the ethanol number density, and (d) water−water, water−ethanol, and water−(ethanol + water), denoted as OW−HW, OW−HE, and OW−HX, plotted against the water number density. The ethanol concentration of the mixtures is given on each figure to aid identification. The standard deviations are also given for each point.

Figure 4 presents calculated coordination numbers as the function of the number density (number density values can be found in Table 1, so that the connection with a given calculation may be made) but the composition of the mixtures is indicated in the figures as well. The standard deviations of the coordination numbers are also represented as error bars in Figure 4 for each data point but due to their small value (maximum is 0.097) they are hardly visible. In Figure 4a, the first data points at ethanol number density ρE = 0.0047 Å−3 correspond to simulations with a 20% ethanol content, the last points at ρE = 0.0103 Å−3 represent pure ethanol; in panel b the first points at water number density ρW = 0.0025 Å−3 correspond to 20% water, and the last points at ρW = 0.0309 Å−3 represent pure water. Interestingly, varying the water force field caused visible (although not excessive) deviations for mixtures with higher water contents not only in partial coordination numbers involving water but in the purely ethanol-related coordination number, as well. The other notable feature is that the coordination number does not change linearly with the number density of the central particles (neither with the density of the neighbors, not shown). In Figure 4a, the total coordination number of oxygen around ethanol oxygen (OE−OX) (X denoting both ethanol and water) changes in a fashion close to linear from 2.01 ± 0.01 for pure ethanol to 2.6−2.8 for Et20, depending on the force field. For pure ethanol, Benmore et al.46 estimated a value of 2 for the OE−OE coordination number from a combined neutron diffraction and MD study, while Noskov et al.20 found 1.94 from their MD study. From the curvatures of the OE−OE and OE−OW graphs it seems that ethanol prefers water to ethanol as first neighbor (similarly to the finding of Noskov et al.20). In Figure 4b, similar trends can be found for the first oxygen neighbors of water oxygen atoms. The total first coordination

HB conditions, and for EHD (STRICT) the TIP4P-2005 force field produced the smallest value. Overall, summing the χ2 for all the concentrations and dimer types for a given force field, TIP4P-2005 produced the closest agreement in case of both the MP4 and the CCSD(T) ab initio data for both the LOOSE and STRICT H-bond definitions (see Table 5). These results strengthen our earlier finding, based on the comparison between experimental and MD calculated X-ray structure factors, that for Et80 the SWM4-DP water force field and for the Et60-Et20 range the TIP4P-2005 water force field give the best agreement. In the case of the 60% ethanol concentration, SWM4-DP provides an equally good description of the system based on the structure factor analysis; this is not supported by the binding energy calculations. 3.4. Coordination Numbers. The ethanol−ethanol, ethanol−water, water−water, and water−ethanol (the first molecule being in the central one) O−O first coordination numbers (ncoord) have been determined from the appropriate partial g(r) functions, as shown in Figure 4. The integration was carried out to the first minimum of the appropriate (intermolecular) g(r), which varied between 3.6 (pure ethanol) to 3.25 Å (Et20 and Et40 with SWM4-DP). Simulations with higher ethanol contents produced slightly larger first coordination spheres. In general, the SPC/E water potential brought about the bulkiest, whereas the SWM4-DP force field produced the smallest coordination spheres among calculations for the same concentration. (See Supporting Information Table 1 for the coordination number ranges.) Interestingly, all the pure water simulations have 3.4 Å as the upper boundary for the first O−O coordination sphere, and their coordination numbers are identical at ∼4.6. That is, the presence of ethanol influences differently the coordination spheres for the different water force fields. G

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. Average number of different kinds of H-bonds/molecule with standard deviations for the (a) TIP4P-2005 and (b) SWM4-DP MD simulations. X denotes either ethanol or water, that is, OE−HE---OX denotes all types of H-bonded pairs containing ethanol hydroxyl hydrogen; E---X means the average number of all kinds of H-bonds/ethanol molecules, while W---X is the average number of all kinds of H-bond/water molecules.

case of the TIP4P-2005 potential, the “participation ratio” has slightly decreased with increasing water concentration (OE− HE---OX, OW−HW---OX) (X denoting either ethanol or water), while for the SWM4-DP force field it has slightly increased, reaching 0.99 for Et20. (The SPC/E force field, not shown, displayed similar behavior as the TIP4P-2005 in tendency, but to a somewhat lesser extent.) The number of H-bonds/ethanol molecule formed through the ethanol oxygen (OE---HX−OX) was 1.1 for the Et80 simulations and increased to 1.4−1.6 for the Et20 simulations (see Figure 5). The number of H-bonds per water molecules formed through water oxygen (OW---HX−OX) was ∼1.4 for Et80 and increased to 1.7−1.8 for the Et20 simulations. This means that in the Et80 simulations the majority of the oxygen nonbonding electron pairs involved in H-bonding belong to ethanol molecules (cf. the much larger number of ethanol molecules, 1226 EtOH versus 316 H2O). Also, only a small ratio of the oxygen atoms (either belonging to ethanol or water) had both nonbonding electron pairs involved in Hbonding. On the other hand, 40% of the water oxygens used both their nonbonding electron pairs. Similarly, in the case of the 20% ethanol containing samples more water molecules were involved in H-bonding through both nonbonding electron pairs than ethanol. The reason for that might be that the nonbonding pairs of the ethanol oxygen are sterically more hindered, due to the presence of the ethyl group. The total average number of H-bonds/water molecule was ∼3.3 for the Et80 simulations, which increased to 3.6−3.8 for Et20, depending on the force field; the corresponding value was ∼3.9 for all the three pure water simulations. The total average number of H-bonds/ethanol molecules changed from 1.9 ± 0.01 for pure ethanol to 2.3−2.6 for Et20, depending on the force field. The number of H-bonds found in a calculation depends on many factors, like the simulation method, the H-bond definition, the applied cutoffs for the O---O and O---H distances and the O−H---O cutoff angle. Narten et al.42 estimated 1.8 H-bonds/ethanol molecule by least-squares fitting to the large-k portion of the k-weighted distinct structure functions. Noskov20 found 3.03 H-bonds/molecule for pure water and 1.65 for pure ethanol in their study of ethanol−water mixtures (using a polarizable ethanol potential combined with the SWM4-DP water force field), applying an O---H cutoff of 2.4 Å and an O−H---O angle >150° as criteria for H-bonding. Oleinikova et al.47 obtained 3.26 H-bonds per molecule in their pure water MD simulation with the TIP4P force field using 3.5 Å for the O---O cutoff and an energy restriction of Uint < −2.6

number curves have different curvatures while changing from ∼3.5 ± 0.05 for models containing 20% water to 4.62−4.64 ± 0.02 in case of pure water for the different water force field simulations. Values for the 80% water models using TIP4P2005 and SWM4-DP water potentials are lower than the one for SPC/E; this originates from the different water−water coordination numbers, as the water−ethanol ncoord is virtually identical for each simulation at a given concentration. The “cumulative” OW−OX curves are not far from linear. Water− water coordination numbers for simulations with higher water content show that water molecules prefer to have more water than ethanol in their first coordination sphere, similarly to ethanol (which agrees, again, with the findings of Noskov et al.20). This is not surprising, as in addition to the higher probability of finding a water molecule anywhere in the simulation box, water is a lot smaller than ethanol, and can fit better to the first coordination sphere. Comparing these findings to the oxygen−hydrogen first coordination numbers (see the Supporting Information Table 2 for distance ranges), the curves in panels a and c of Figure 4 are similar to those in panels b and d, respectively. It has to be noted that, especially for hydrogen around ethanol oxygen, there is hardly any difference between the curves belonging to different force fields at the same concentration. The OW−HW coordination numbers and, consequently, the total OW−HX ones (X denoting either water or ethanol) do not differ as much as the corresponding OW−OW values: even the tendency is different, SPC/E ncoord (OW−HW) being the lowest out of the simulations with 80% water while it possesses the largest ncoord (OW−OW). 3.5. Number of Hydrogen Bonds. The number of hydrogen bonds for both the LOOSE and STRICT conditions (H-bond definition are given in Section 2.2) were evaluated for the different components separately; some of the results are shown in Figure 5 for simulations using the TIP4P-2005 and SWM4-DP force fields. The number of H-bonds for the STRICT conditions is almost identical (∼99%) to the values calculated for the LOOSE conditions for every concentration and type, so only values obtained using the LOOSE condition will be discussed. The total number of H-bonds increases with increasing water concentration, as expected (not shown); a water molecule is theoretically capable of forming four H-bonds while an ethanol molecule only three. Among the three force fields, the SWM4DP force field produces the highest total number of H-bonds, although differences are not significant. It is interesting to notice that for all concentrations and force fields, almost all hydroxyl hydrogen atoms were involved in H-bonding. In the H

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

produced the strongest first maximum of the Wat−Wat CMD. At medium and high water concentrations, the polarizable SWM4-DP force field clearly produces sharper distributions than the other two force fields do. Note that results from this potential agree the least with experimental X-ray data over the medium and high water concentration ranges and therefore, the sharpness of the corresponding CMDs cannot be considered as a genuine feature of the structure. Regarding the second maximum (see insert in Figure 6), both the peak heights and positions change with concentration. Considering results belonging to the simulations possessing the best agreement with the experimental data at each concentration, the second peak position has shifted from 5.05 Å (Et80 SWM4-DP) to 4.85 Å (Et60 SWM4-DP, TIP4P-2005), 4.75 Å (Et40 TIP4P-2005), and 4.65 Å (Et20 TIP4P-2005). That is, increasing the water concentration does not only increase the number of second neighbors, but at the same time the second coordination sphere becomes more compact. There are more significant differences between the ethanol− water CMDs, both concerning distributions at the same concentrations and different force fields and CMDs with the same force fields at different concentrations (Figure 7).

kcal/mol. Soper et al.48 calculated 3.58 H-bonds/water molecule using an O---O cutoff value of 3.5 Å and O−H---O angle >150° in their work based on analyses of their neutron diffraction data. In our present work, cutoff values for water O---O and O---H of ∼3.4 and 2.5 Å, respectively, were applied (in addition with an angular restriction of O−H---O > 120° for the STRICT condition based on the definition of Chen et al.7 and on our previous work33). In order to check how sensitive the value of the number of H-bond/molecule is to the angular restriction, the number of H-bonds/molecule was also calculated using O−H---O > 150° for the three water models considered here. Under these conditions, the number of Hbonds/water molecule was 3.0 ± 0.01 (SPC/E) and 3.1 ± 0.01 (TIP4P-2005 and SWM4-DP), which values are very similar to results of Noskov et al.20 3.6. Center-of-molecule Distribution. The “center-ofmolecule distributions” (CMD) were calculated for all the possible combinations, Et−Et, Et−Wat, and Wat−Wat. The“center-of-molecule” is defined as the geometric center of a given molecule; that is, it is essentially the center-of-mass, COM, assuming a uniform mass for each atom. This construction was chosen because the center-of-molecule represents better the space occupation of a molecule than the center-of-mass would be able to. The ethanol−ethanol CMDs are relatively similar, regardless of the force field and concentration (not shown); notable differences are that the position of the first maximum for all the three force fields shifts from 4.65 to 4.85 Å with decreasing ethanol concentration, while the second peak moved to smaller distances from 9.05 to 8.55 Å. These changes show that the presence of water does affect the first and second Et−Et coordination shells and, moreover, that the changes are not due to pure steric effects. Concerning the water−water CMD (Figure 6), the position of the first maximum is 2.85 Å for all the force fields at each

Figure 7. Ethanol−water CMDs for the SPC/E, TIP4P-2005 (shifted along the y-axis by 1), and SWM4-DP (shifted by 2) water potentials at 20, 40, 60, and 80% ethanol concentration (in the order of increasing peak heights).

Inspecting results from the SWM4-DP force field, there are two distinct peaks in the region of 2.75−5.35 Å at 3.25 and 4.15 Å at each concentration; the degree of splitting increases with the decreasing ethanol concentration. For the TIP4P-2005 and the SPC/E water models, the two maxima are merged and the distinction between the two “submaxima” becomes harder with the decreasing ethanol concentration. For analyzing the peak structure further, center-of-molecule distributions have been calculated separately for H-bonded and non-H-bonded pairs of molecules (LOOSE condition; definitions for H-bonding are given in Section 2.2). For the mixed, ethanol−water, distribution two possibilities are distinguished, depending on whether the H atom of the Hbond belongs to an ethanol or a water molecule. Interestingly, there is a significant difference between CMDs of the OE−HE--O W and the O W −H W ---O E H-bonded molecule pairs, independently of the concentration and the force field applied. As an example, the total Et−Wat, as well as the H-bonded and not H-bonded CMDs are shown in Figure 8 for the Et80 SWM4-DP and Et20 TIP4P-2005 simulations, which provided the best agreement with experimental data for the given concentrations. If the H atom in a mixed (Et−Wat) H-bond belongs to the water molecule, then the center-of-molecule distribution has only one peak around 3.25 Å; a characteristic conformation that contributes to this maximum is given in Figure 8. This

Figure 6. Water−water center-of-molecule distributions for all MD simulations; the region of the second maximum is enlarged in the inset. The simulations providing the best fit to experiment at a given concentration are emphasized by thick lines and bold characters.

concentration, but there are significant differences between the force fields in terms of the heights of the first and second peaks belonging to the same concentration. The TIP4P-2005 force field produces the lowest first maximum at each concentration. The first peaks produced by the SWM4-DP force field are the highest for the 20−60% ethanol range, while at the highest (80%) ethanol concentration, it is the SPC/E force field that I

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

ethanol molecule is bent toward the water molecule are the most frequent. The relation between the Et−Wat H-bonded pair distributions is nearly uniform at each concentration; the different peak shapes are results of differences in terms of the distributions for the non-H-bonded pairs. These pairs consist of ethanol and water molecules that, although close enough to each other, are in the wrong orientation to form a hydrogen bond. In the ethanol-rich systems, there are more of these pairs with CMD ∼3.25 Å than in the water-rich systems. 3.7. Conformational Analyses Based on Bivariate (Distance−Angle) Distributions. Two dimensional distributions C[r, cos(θ)], depending on the distance between the two molecular centers and the cosine of the O−H---O H-bond angle cos(θ), for H-bonded molecule pairs have also been calculated. The relationship between C[r, cos(θ)] and the CMD distribution is the following: CMD(r) = Σθ180 =0 C[r, cos(θ)]. These kinds of distributions have not revealed any interesting inner structure (therefore they are not shown); the one maximum in each decays quickly along the cos(θ) axis. On the other hand, it is worth quoting some statistics. About 29% of the H-bonded molecule pairs (applying the LOOSE condition) had H-bond angles between 169 and 180° for all combinations (OE−HE---OE, OE−HE---OW, OW−HW---OE, OW−HW---OW) in mixture with 80% ethanol in simulations with the SPC/E and TIP4P-2005 force fields; this ratio decreased to ∼24% for samples with 20% ethanol. During this transition, the maximum has become a little more diffuse: 99% of the pairs had H-bond angle >120° in Et80, which decreased to ∼98% in Et20 simulations (see Figure 9a; only results for the TIP4P-2005 water force field are shown). That is, the SPC/E and TIP4P-2005 force fields influenced similarly not only the Wat−Wat and Wat−Et, but the Et−Et combination, as well. Interestingly, calculations with the SWM4-DP potential provide a markedly different outcome (see Figure 9b): the C[r, cos(θ)] distribution becomes sharper with the decreasing ethanol concentration for all combinations (even for the pure Et−Et one). The effect increases with the involvement of water in the H-bonds in the order OE−HE---OE < OE−HE---OW
169° for (a) the TIP4P-2005 and (b) the SWM4-DP water force fields. J

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 10. Bivariate distribution C[rCB‑CB, cos(θ)], depending on the distance of the methyl C atoms of the H-bonded ethanol pairs and the H-bond angle, for the Et80 SWM4-DP simulation. The conformation of the ethanol pairs characteristic to the peak regions are also shown. H atoms are light gray, C atoms are dark gray, O atoms are red, and the H-bond is depicted by yellow lines.

OW−HW---OE < OW−HW---OW. Concerning the SWM4-DP Et20 simulation, 42% of the H-bonded molecules of the Et−Et and 70% of the Wat−Wat pairs had rather straight, >169° Hbond angle; on the other hand, the TIP4P-2005 Et20 calculation yields only 27% of the H-bonded pairs of molecules had a >169° H-bond. At the highest ethanol concentration, in the SWM4-DP Et80 simulation, 99% of the molecular pairs had H-bonds >120° H-bond for each combination and even this value has increased to ∼100% for the Et20 case. Clearly, the SWM4-DP water force field due to its polarizable nature has a significant influence on the H-bond angle, resulting in straighter H-bonds, even in case of the Et−Et combination. This feature may prove to be useful for mixtures with higher ethanol content, as for Et80 the SWM4-DP simulation produced the overall best agreement with the X-ray data, and for Et60 the SWM4-DP and TIP4P-2005 simulations produced equally good agreement. However, in mixtures with less ethanol the polarizable water potential seems to overestimate the regularity (i.e., the “straightness”) of hydrogen bonds. In order to further analyze the conformations of the Hbonded Et−Et pairs, two-dimensional distributions C[rCB‑CB, cos(θ)], depending on the distance between the methyl C atoms (rCB‑CB) and the cosine of the O−H---O angle, cos(θ), have been calculated. The relationship between C[rCB‑CB, cos(θ)] and the PRDF for the methyl C atoms, g(rCB‑CB), is the following: g(rCB−CB) = Σθ180 =0 C[rCB−CB, cos(θ)], similarly to the center-of-molecule distribution. For all the simulations, regardless of the composition and the water force field (if applied), the distribution has two main peaks, located at θ = 180° and rCB‑CB = ∼4.2 and 5.2 Å, as can be seen in Figure 10 (shown for the Et80 SWM4-DP simulation, as an example). The conformations of several molecule pairs were visually inspected. This process revealed that (a) the first, less intense, peak around 4.2 Å is formed by ethanol pairs where both the methyl groups are inclined toward the other ethanol molecule or perpendicular to the H-bond, and that (b) in the majority of H-bonded Et−Et pairs one of the methyl groups is inclined toward the other molecule or perpendicular to the H-bond, while the other one points away from the molecular center of

the pair. (For instructive drawings of these conformations, see also Figure 10.) With increasing water concentration the maxima connected to the SWM4-DP simulations become narrower and more intense than the peaks from the other two water models at the same concentration. This effect is excessive: for instance, for the Et 20 SWM4-DP simulation peak heights are three times of what result from the SPC/E and TIP4P-2005 simulations. As it was discussed earlier in connection with Figure 5, the number of ethanol−ethanol H-bonds is nearly equal for the SPC/E and SWM4-DP force fields and only 15% smaller for the TIP4P2005 potential, so this cannot account for the difference. It has to be attributed to the polarizability of the SWM4-DP force field; it is interesting that the choice of the water force field can have such a large effect on the ethanol−ethanol distribution, as well. Similarly, distributions C[rCB‑O, cos(θ)] of the H-bonded pairs, depending on the distance between the methyl C atom (CB) and the oxygen atom of the other molecule and the cosine of the O−H---O H-bond angle, have been calculated. The relationship between C[rCB‑O, cos(θ)] and the PRDF for the methyl C atom−O atom, g(rCB‑O), is the following: g(rCB−O) = 180 C[rCB−O, cos(θ)], similarly as before. Four different Σθ=0 distributions were considered separately, depending on the type of the molecule and whether the CB atom belongs to the hydrogen donor in the hydrogen bond. The variables of the four distributions are the following: the distance between the methyl C atom, CB, of an H-donor Et molecule and the O atom of the H-bonded molecule, plus the H-bonding angle, cos(θ) (denoted as CBHD) for (1) ethanol−ethanol [H3CB−CH2− O−H---O(H)−CH2CH3] and (2) ethanol−water [H3CB− CH2−O-H---OH2] pairs; the distance between the methyl C atom, CB, of an H-acceptor Et molecule and the O atom of the H-bonded molecule, plus the H-bonding angle, cos(θ) (denoted as CBHA) for (3) ethanol−ethanol [H3CB−CH2− (H)O---H−O−CH2CH3] and (4) ethanol−water [H3CB− CH2−(H)O---H−O−H] pairs (the bold, italic characters denote the atoms involved in the distance distribution calculation). K

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 11. Distributions of the distances between the methyl C atom (CB) of an ethanol molecule and the O atom of another ethanol molecule, g(rCB‑O), for the Et80 SWM4-DP simulation. CB belongs (a) to an H-donor ethanol (CBHD) and (b) to an H-acceptor ethanol molecule (CBHA). (c) Bivariate C[rCB‑O, cos(θ)] distributions for the CBHD-type H-bonded ethanol pairs; (d) similarly for the CBHA pairs. Atom and H-bond coloring is the same as in Figure 10. Characteristic conformations are also shown; dashed lines indicate the CB−O distances in question.

Figure 12. Distributions of the distances between the methyl C atom (CB) of an ethanol molecule and the O atom of a water molecule, g(rCB‑O), for the Et80 SWM4-DP simulation. CB belongs (a) to an H-donor ethanol (CBHD) and (b) to an H-acceptor ethanol molecule (CBHA).

depending on the CB−O distance can be detected, apart from the straight angle. Characteristic conformations that contribute to the peak maxima are also shown. With the decreasing ethanol concentration, the separation of the two peaks in the CBHD-type distribution is increasing. Peaks are higher and narrower for the SWM4-DP force field than for the other two potentials, similarly to the CB−CB distributions. In case of the ethanol−water g(rCB‑O) pairs (Figure 12), distributions for the H-bonded pairs are similar to that of the corresponding ethanol−ethanol distributions. Interestingly, the total distributions are different; while for the ethanol−water

First looking at the g(rCB-O) distributions for the H-bonded ethanol−ethanol molecule pairs, it is interesting to notice that while for the CBHD-type H-bonded pairs there are two peaks at ∼3.7 and 4.8 Å, both contributing to the cumulative distribution (Figure 11a), the distribution for the CBHA-type H-bonded ethanol−ethanol pairs is different, showing only one maximum at 3.7 Å (and only a small shoulder around 4.8 Å, Figure 11b). That is, in this case the peak at 4.8 Å on the total distribution can mainly be attributed to ethanol pairs not involved in this type of H-bonding. The corresponding C[rCB‑O, cos(θ)] distributions for the H-bonded pairs are shown in Figure 11c,d: no preference for any specific H-bond angle L

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B pairs the first peak at ∼3.7 Å is larger and the second one at ∼4.8 Å is smaller, for the ethanol−ethanol pairs the first maximum is less intense than the second one. This can be explained by the fact that ethanol−water pairs can get closer to each other more easily, due to the smaller size of the water molecule. Bivariate C[rCB‑O, cos(θ)] distributions were also calculated as before; they are very similar to the corresponding ethanol− ethanol distributions (and therefore, not shown). Characteristic conformations for H-bond angles ∼180° are also given in Figure 12. Features of these distributions hardly change with decreasing ethanol concentration and, similarly to the previous cases, calculations with the SWM4-DP water potential produce sharper and higher maxima than the two other force fields.

the simulations performed. Water oxygen atoms are more often participate in hydrogen bonding using both of their nonbonding electron pairs than ethanol oxygen atoms, most probably due steric considerations. The average number of H-bonds/water molecule for the mixture with 80 mol % ethanol was ∼3.2, which increased to 3.6−3.8 for the 20 mol % ethanol solution, depending on the force field; the value for all the three pure water simulations was ∼3.8. The average number of H-bonds/ethanol molecules was 1.9 for pure ethanol, which increased to 2.3−2.6 for the 20 mol % ethanol systems, depending on the force field. The OW−OW first coordination number for pure water was ∼4.6, which is significantly higher than that of ∼3.8 H-bond/ water molecule. This shows that the assumption that the number of H-bonds may satisfactorily be estimated from the coordination number41 is questionable. The same conclusion can be drawn from the comparison of the other (Et−Et, Et− Wat) types of H-bonds to the appropriate O−O coordination number. It is worth emphasizing again that the number of hydrogen bonds per molecule strongly depends on several factors, such as the calculation method and the hydrogen bond definition. Comparing values coming from different studies with different hydrogen bond definitions is therefore not advisable: systematic series of values (coming possibly from the same study) may better reveal important trends. Our work suggests that the basis of the hydrogen bond definition should rather be the “two distance criterium” (LOOSE condition), and the use of the additional (arbitrary) 120° angular constraint (STRICT condition) is superfluous. Center-of-molecule distributions revealed that the ethanol− ethanol first coordination sphere expands from 4.65 to 4.85 Å with increasing water concentration, whereas the second sphere contracts from 9.05 to 8.55 Å. The size of the first water−water coordination sphere does not change, while the second sphere also contracts from 5.05 to 4.65 Å. This is an indication that both the ethanol and water H-bonded networks become denser with the increasing water concentration. The mixed ethanol− water CMD revealed that if the H-bond is formed through the water hydrogen, then the centers of ethanol and water molecule are located closest to each other, around ∼3.25 Å; on the other hand, if the H-bond is formed through the ethanol hydroxyl hydrogen then centers are found further away, around 3.65 and 4.25 Å. The H-bond angle is strongly influenced by the force field even for the ethanol−ethanol pairs. The SWM4-DP water force field systematically provided much sharper angular distributions than the other two potentials. As the SWM4-DP fails to yield a good description of the experimental data at lower ethanol concentrations, it seems that the high ratio of nearly straight hydrogen bonds is not a characteristic feature of the real system. The ratio of ∼27% suggested by the TIP4P-2005 force field may be a better guess, as this potential provides the best agreement with experimental X-ray data. The distribution of the distance between the methyl carbon atom and the O atom of its hydrogen-bonded pair differs for the following two cases: (1) the methyl carbon atom belongs to the hydrogen donor (CBHD), and (b) the methyl carbon atom belongs to the hydrogen acceptor ethanol molecule (CBHA). So far only orientations of the hydrogen-bonded molecule pairs have been analyzed; work is in progress in order to investigate the hydrogen-bonded network in detail, including

4. CONCLUSIONS A series of molecular dynamics simulations has been performed for 20, 40, 60, and 80 mol % ethanol-in-water mixtures, pure ethanol, and pure water using the OPLS-AA force field for ethanol and the SPC/E, TIP4P-2005, and SWM4-DP potential models for water. The total scattering X-ray structure factor was applied for assessing the applicability of these interaction potentials. No single water force field could be found to provide the best agreement between the experiment and simulation for all concentrations. For the mixture with 80% ethanol the SWM4-DP, with 60% ethanol the SWM4-DP and TIP4P-2005, whereas at lower ethanol concentration the TIP4P-2005 water force field proved to be the most successful. The binding energies of the hydrogen-bonded ethanol− water pairs were also calculated for both the LOOSE and STRICT HB conditions and compared to gas-phase MP4 and CCSD(T) ab initio results of Fileti.27 A generally good agreement could be found, regardless of the fact that hydrogenbonded pairs in the liquid phase not necessarily have the same optimum geometry as in gas phase, not to mention that in the present mixtures the hydrogen-bonded pairs are not isolated, but mostly parts of larger clusters. Overall, the binding energy results support the findings based on the comparison of the calculated and experimental structure factors. Partial radial distribution functions have been calculated and examined in order to reveal major differences between potential models and also as a function of composition. The second maximum of the water−water partial radial distribution function that is supposed to be the signature of the hydrogen-bonded network41 shifts toward higher distances with increasing ethanol concentration; this is in contradiction with the finding of Dixit et al.,41 for a 70 mol % methanol− water mixture. This indicates that the larger ethanol molecule disrupts the water structure significantly, as compared to the smaller methanol molecules. The number of oxygen atoms around ethanol oxygen (OE− OX) grows in a fashion close to linear from ∼2 (pure ethanol) to 2.6−2.8 (20 mol % ethanol), depending on the force field. The curvatures of the coordination number versus number density curves indicate that ethanol molecules prefer water to ethanol as first neighbor. Having investigated oxygen first neighbors around water oxygen atoms it could be established that water molecules also prefer to have more water than ethanol in their first coordination spheres. The number of H-bonds/molecules for all the possible molecule pair combination was determined. Almost all the hydroxyl hydrogen atoms participate in hydrogen bonding in all M

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(12) Wensink, E. J. W.; Hoffmann, A. C.; van Maaren, P. J.; van der Spoel, D. Dynamic Properties of Water/Alcohol Mixtures Studied by Computer Simulation. J. Chem. Phys. 2003, 119, 7308−7317 DOI: 10.1063/1.1607918. (13) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. J. A New Force Field for Molecular Mechanical Simulation of Nucleic Acids and Proteins. J. Am. Chem. Soc. 1984, 106, 765−784 DOI: 10.1021/ja00315a051. (14) Jorgensen, W. L. Optimized Intermolecular Potential Functions for Liquid Alcohols. J. Phys. Chem. 1986, 90, 1276−1284 DOI: 10.1021/j100398a015. (15) Van Leeuwen, M. E. Prediction of the Vapour-Liquid Coexistence Curve of Alkanols by Molecular Simulation. Mol. Phys. 1996, 87, 87−101 DOI: 10.1080/00268979600100031. (16) Gotlib, I.Yu.; Piotrovskaya, E. M. Properties of Coexisting Phases for the Ethanol-Ethane Binary System by Computer Simulation. J. Phys. Chem. B 1999, 103, 7681−7686 DOI: 10.1021/ jp990263j. (17) Kaminski, G.; Duffy, E. M.; Matsui, T.; Jorgensen, W. L. Free Energies of Hydration and Pure Liquid Properties of Hydrocarbons from the OPLS All-Atom Model. J. Phys. Chem. 1994, 98, 13077− 13081 DOI: 10.1021/j100100a043. (18) Jorgensen, W. L.; Maxwell, D.; Tirado-Rives, S. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236 DOI: 10.1021/ja9621760. (19) Gao, J.; Habibollazadeh, D.; Shao, L. A Polarizable Intermolecular Potential Function for Simulation of Liquid Alcohols. J. Phys. Chem. 1995, 99, 16460−16467 DOI: 10.1021/j100044a039. (20) Noskov, S. Yu.; Lamoureux, G.; Roux, B. Molecular Dynamics Study of Hydration in Ethanol-Water Mixtures Using a Polarizable Force Field. J. Phys. Chem. B 2005, 109, 6705−6713 DOI: 10.1021/ jp045438q. (21) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271 DOI: 10.1021/j100308a038. (22) Tarek, M.; Tobias, D. J.; Klein, M. L. Molecular Dynamics Investigation of an Ethanol-Water Solution. Physica A 1996, 231, 117− 122 DOI: 10.1016/0378-4371(95)00450-5. (23) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935 DOI: 10.1063/1.445869. (24) Lamoureux, G.; MacKerell, A. D., Jr.; Roux, B. A Simple Polarizable Model of Water based on Classical Drude Oscillators. J. Chem. Phys. 2003, 119, 5185−5197 DOI: 10.1063/1.1598191. (25) Bernal, J. D.; Fowler, R. H. A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl Ions. J. Chem. Phys. 1933, 1, 515−548 DOI: 10.1063/1.1749327. (26) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505/1−12 DOI: 10.1063/1.2121687. (27) Fileti, E. E.; Chaudhuri, P.; Canuto, S. Relative Strength of Hydrogen Bond Interaction in Alcohol−Water Complexes. Chem. Phys. Lett. 2004, 400, 494−499 DOI: 10.1016/j.cplett.2004.10.149. (28) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718 , http://www.gromacs.org DOI: 10.1002/jcc.20291. (29) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092 DOI: 10.1063/1.464397. (30) Essman, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593 DOI: 10.1063/1.470117. (31) (a) Nose, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268 DOI: 10.1080/00268978400101201. (b) Hoover, W. G. Canonical

the occurrence of ring structures in ethanol−water systems of various compositions.

■

ASSOCIATED CONTENT

S Supporting Information *

The oxygen−oxygen and oxygen−hydrogen nonbonding distance ranges for the hydrogen bond determination are provided for all the simulations. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +36 30 273 5516. Fax: +36 1 392 2215. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. László Temleitner (Wigner RCP, Budapest, Hungary) for providing his unpublished X-ray diffraction data with us. The data were taken at the SPring-8 Synchrotron Facility (Hyogo, Japan) under beamtime allocations 2011B1554 and 2013A1083. This work was supported by the National Basic Research Fund (OTKA; Hungary), under contract No. 083529.

■

REFERENCES

(1) Stewart, G. W.; Morrow, R. M. X-Ray Diffraction in Liquids: Primary Normal Alcohols. Phys. Rev. 1927, 30, 232−244 DOI: 10.1103/PhysRev.30.232. (2) Stewart, G. W.; Skinner, E. W. X-Ray Diffraction in Liquids: A Comparison of Certain Primary Normal Alcohols and their Isomers. Phys. Rev. 1928, 31, 1−9 DOI: 10.1103/PhysRev.31.1. (3) Raman, C. V.; Sogani, C. M. X-Ray Diffraction in Liquids. Nature 1927, 119, 601−601 DOI: 10.1038/120514a0. (4) Zhang, C.; Yang, X. Molecular Dynamics Simulation of Ethanol/ Water Mixtures for Structure and Diffusion Properties. Fluid Phase Equilib. 2005, 231, 1−10 DOI: 10.1016/j.fluid.2005.03.018. (5) Laaksonen, A.; Kusalik, P. G.; Svishchev, I. M. ThreeDimensional Structure in Methanol-Water Mixtures. J. Phys. Chem. A 1997, 101, 5910−5918 DOI: 10.1021/jp970673c. (6) Zhang, N.; Li, W.; Chen, C.; Zuo, J.; Weng, L. Molecular Dynamics Study on Water Self-Diffusion in Aqueous Mixtures of Methanol, Ethylene Glycol and Glycerol: Investigations from the Point of View of Hydrogen Bonding. Mol. Phys. 2013, 11, 939−949 DOI: 10.1080/00268976.2012.760050. (7) Chen, B.; Siepmann, J. I. Microscopic Structure and Solvation in Dry and Wet Octanol. J. Phys. Chem. B 2006, 110, 3555−3563 DOI: 10.1021/jp0548164. (8) Vrhovšek, A.; Gereben, O.; Pothoczki, Sz.; Tomšič, M.; Jamnik, A.; Kohara, S.; Pusztai, L. An Approach towards Understanding the Structure of Complex Molecular Systems: the Case of Lower Aliphatic Alcohols. J. Phys.: Condens. Matter 2010, 22, 404214/1−9 DOI: 10.1088/0953-8984/22/40/404214. (9) Harris, K. R.; Newitt, P. J. Diffusion and Structure in Dilute Aqueous Alcohol Solutions: Evidence for the Effects of Large Apolar Solutes on Water. J. Phys. Chem. B 1998, 102, 8874−8879 DOI: 10.1021/jp9820370. (10) Petong, P.; Pottel, R.; Kaatze, U. Water-Ethanol Mixtures at Different Compositions and Temperatures. A Dieletric Relaxation Study. J. Phys. Chem. A 2000, 104, 7420−7428 DOI: 10.1021/ jp001393r. (11) Tan, M.-L.; Cendagorta, J. R.; Ichiye, T. Effects of Microcomplexity on Hydrophobic Hydration in Amphiphiles. J. Am. Chem. Soc. 2013, 135, 4918−4921 DOI: 10.1021/ja312504q. N

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697 DOI: 10.1103/PhysRevA.31.1695. (32) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962 DOI: 10.1002/jcc.540130805. (33) Vrhovšek, A.; Gereben, O.; Jamnik, A.; Pusztai, L. Hydrogen Bonding and Molecular Aggregates in Liquid Methanol, Ethanol, and 1-Propanol. J. Phys. Chem. B 2011, 115, 13473−13488 DOI: 10.1021/ jp206665w. (34) Temleitner, L. Personal communication, 2013. (35) Kohara, S.; Suzuya, K.; Kashihara, Y.; Matsumoto, N.; Umesaki, N.; Sakai, I. A Horizontal Two-Axis Diffractometer for High-Energy XRay Diffraction using Synchrotron Radiation on Bending Magnet Beamline BL04B2 at SPring-8. Nucl. Instrum. Methods Phys.Res. A 2001, 467/468, 1030−1033 DOI: 10.1016/S0168-9002(01)00630-1. (36) Temleitner, L. Structure Determination of Liquid Carbon Tetrabromide via a Combination of X-Ray and Neutron Diffraction Data and Reverse Monte Carlo Modeling. J. Mol. Liq. 2014, 197, 204− 210 DOI: 10.1016/j.molliq.2014.05.011. (37) Kohara, S.; Itou, M.; Suzuya, K.; Inamura, Y.; Sakurai, Y.; Ohishi, Y.; Takata, M. Structural Studies of Disordered Materials using High-Energy X-Ray Diffraction from Ambient to Extreme Conditions. J. Phys.: Condens. Matter 2007, 19, 506101/1−15 DOI: 10.1088/09538984/19/50/506101. (38) Hart, R. T.; Benmore, C. J.; Neuefeind, J.; Kohara, S.; Tomberli, B.; Egelstaff, P. A. Temperature Dependence of Isotopic Quantum Effects in Water. Phys. Rev. Lett. 2005, 94, 047801/1−4 DOI: 10.1103/ PhysRevLett.94.047801. (39) Pusztai, L.; Pizio, O.; Sokolowski, S. Comparison of Interaction Potentials of Liquid Water with respect to their Consistency with Neutron Diffraction Data of Pure Heavy Water. J. Chem. Phys. 2008, 129, 184103/1−6 DOI: 10.1063/1.2976578. (40) Steinczinger, Zs; Pusztai, L. Comparison of the TIP4P-2005, SWM4-DP and BK3 Interaction Potentials of Liquid Water with Respect to their Consistency with Neutron and X-Ray Diffraction Data of Pure Water. Cond. Matt. Phys. 2013, 16, 43604/1−8 DOI: 10.1063/ 1.2976578. (41) Dixit, S.; Crain, J.; Poon, W. C. K.; Finnley, J. L.; Soper, A. K. Molecular Segregation Observed in a Concentrated Alcohol-Water Solution. Nature 2002, 416, 829−832 DOI: 10.1038/416829a. (42) Narten, A. H.; Habenschuss, A. Hydrogen Bonding in Liquid Methanol and Ethanol Determined by X-ray Diffraction. J. Chem. Phys. 1984, 80, 3387−3391 DOI: 10.1063/1.447093. (43) Kakar, R. K.; Quade, C. R. Microwave Rotational Spectrum and Internal Rotation in Gauche Ethyl Alcohol. J. Chem. Phys. 1980, 72, 4300−4307 DOI: 10.1016/1.439723. (44) van Erp, T. S.; Meijer, E. J. Ab Initio Molecular Dynamics Study of Aqueous Solvation of Ethanol and Ethylene. J. Chem. Phys. 2003, 118, 8831−8840 DOI: 10.1063/1.1567258. (45) Masella, M.; Flament, J. P. Relation Between Cooperative Effects in Cyclic Water, Methanol/Water, and Methanol Trimers and Hydrogen Bonds in Methanol/Water, Ethanol/Water, and Dimethylether/Water Heterodimers. J. Chem. Phys. 1998, 108, 7141−7151 DOI: 10.1063/1.476131. (46) Benmore, C. J.; Loh, Y. L. The Structure of Liquid Ethanol: A Neutron Diffraction and Molecular Dynamics Study. J. Chem. Phys. 2000, 112, 5877−5883 DOI: 10.1063/1.481160. (47) Oleinikova, A.; Brovchenko, I.; Geiger, A.; Guillot, B. Percolation of Water in Aqueous Solution and Liquid−Liquid Immiscibility. J. Chem. Phys. 2002, 117, 3296−3304 DOI: 10.1063/ 1.1493183. (48) Soper, A. K.; Bruni, F.; Ricci, M. A. Site−Site Pair Correlation Functions of Water from 25 to 400°C: Revised Analysis of New and Old Diffraction Data. J. Chem. Phys. 1997, 106, 247−254 DOI: 10.1063/1.473030.

O

DOI: 10.1021/jp510490y J. Phys. Chem. B XXXX, XXX, XXX−XXX