Temperature-dependent structure of methanol-water mixtures on cooling: X-ray and neutron diffraction and molecular dynamics simulations

Methanol-water liquid mixtures have been investigated by high-energy synchrotron X-ray and neutron diffraction at low temperatures. We are thus able to report the first complete sets of both X-ray and neutron weighted total scattering structure factors over the entire composition range (at 12 different methanol concentrations (x$_M$) from 10 to 100 mol%) and at temperatures from ambient down to the freezing points of the mixtures. The new diffraction data may later be used as reference in future theoretical and simulation studies. Measured data are interpreted by molecular dynamics simulations, in which the all atom OPLS/AA force field model for methanol is combined with both the SPC/E and TIP4P/2005 water potentials. Although the TIP4P/2005 water model was found to be somewhat more successful, both combinations provide at least semi-quantitative agreement with measured diffraction data. From the simulated particle configurations, partial radial distribution functions, as well as various distributions of the number of hydrogen bonds have been determined. As a general trend, the average number of hydrogen bonds increases upon cooling. However, the number of hydrogen bonds between methanol molecules slightly decreases with lowering temperatures in the concentration range between ca. 30 and 60 mol % alcohol content. The same is valid for water-water hydrogen bonds above 70 mol % of methanol content, from room temperature down to 193 K.

Methanol, the simplest alcohol, is the closest analogue to water, in which one (hydrophilic) proton has been replaced by a (hydrophobic) methyl group. Despite their similar shape and size, this substitution leads to significant differences in terms of the structure of the two liquids, primarily due to the hydrophobic behavior of the methyl group. While the fourfold coordination of water molecules results in a three-dimensional hydrogen bonded (HB) network in liquid water, methanol molecules participate in at most two hydrogen bonds in practice, and form smaller chains [16][17][18][19]. In aqueous solutions of methanol, both methanol and water can act as H-donor and H-acceptor, and hydrogen bonded networks can be constructed via water-water, methanol-methanol and methanol-water hydrogen bonded pairs.
Most of the publications listed above are limited to ambient temperature: only a much smaller number of studies deal with the temperature dependence of the structure [5,8,20,39,43,52,56,60,61,62,64,65]. The same is valid concerning diffraction: while several studies have considered neutron [41,42,48,49] and X-ray [14,50,51] measurements at ambient temperature, diffraction data at lower temperatures are insufficient. Up to now (and to our knowledge), only one set of X-ray diffraction data is available, for the water-rich solutions (at xM = 0.2, 0.3 and 0.4) [52] only, while two papers are available about neutron diffraction data [5,43] at two methanol concentrations (xM = 0.27 and 0.54), with 2 or 3 temperature points. Interpretation of the X-ray data suggested that the structure of dominant clusters formed in the mixtures at 298 K do not significantly change with lowering the temperatures. The neutron data mentioned above allowed speculations that the 'micro-segregation' between methanol and water clusters is enhanced by lowering the temperature.
The most complete investigation of concentration and temperature dependence of methanol-water solutions was conducted by a series of NMR measurements [61,62,64,65]. Above a concentration dependent crossover temperature hydrophobicity was suggested to play a significant role, while below this temperature the tetrahedral network of water molecules network appeared to determine the properties of the solutions.
In a recent publication molecular dynamics simulation data were compared with X-ray diffraction results of Takamuku et al. [8], where the H-bonding structure have been analyzed from simulated particle configurations. Necessarily, that work was also limited to the water-rich regime. To extend the investigations to equimolar and methanol-rich mixtures would be desirable, but a complete diffraction data set over the whole range of methanol concentration at temperatures down to the freezing point was missing.
By the present study, we provide complementary experimental data on methanol-water liquid mixtures for 12 methanol concentrations, from xM = 0.1 to 1.0. Both X-ray and neutron diffraction total scattering structure factors have been determined at several temperatures, from 300 K down to the freezing points of the mixtures. For a primary interpretation of these, rather large amount of, data,

Experimental details
Both protonated (CH3OH) and deuterated (CD3OD) forms of methyl alcohol (methanol) were purchased from Sigma-Aldrich Co. 12 methanol-water samples were taken for the diffraction experiments, using protonated compounds for X-ray, and the fully deuterated form of methanol and water for neutron measurements. Mixtures with methanol concentrations from 10 mol% to 100 mol%, with 10 mol% equimolar steps, and additionally with 54.42 and 73.37 mol% (peritectic concentration and close to the eutectic concentration, respectively), were prepared. Nominal and exact compositions are shown in Table 1. in the Q-range between 12 and 13 Å -1 , datasets were taken only up to 11.8 Å -1 for further investigations.

X-ray diffraction experiments
Each measurement was controlled by software available at the beamline. The summarized datasets were corrected by efficiency, using a vanadium standard. Intensities from the empty container were also subtracted.
All temperatures and compositions visited by X-ray and neutron diffraction experiments are collected in Table 1, and shown superimposed on the phase-diagram of methanol-water mixtures in Fig. 1. Common features of the XRD structure factors over the entire concentration range are that the first (around Q = 1.85 Å -1 ) and second maxima (around Q = 2.9 Å -1 ) become more intense upon cooling.

Molecular dynamics simulations
At room temperature the height of the 2 nd peak is decreasing with increasing methanol concentration, and only a shoulder can be observed for the xM ≥ 0.8 mixtures. As temperature decreases then a peak around Q = 2.9 Å -1 can be seen even in pure methanol.
Positions of maxima also depend on temperature. At low methanol concentration (xM = 0.1) the position of the first peak moves to smaller Q values as the temperature decreases. As the methanol content increases then the 1st peak position becomes less sensitive to temperature variations --even constant for the composition range 0.5 ≤ xM ≤ 0.7. Increasing methanol concentration further makes the position of the 1 st peak moving to slightly higher Q values as the temperature decreasing. The composition dependence of the position of the 1 st peak is significant at room temperature, but becomes less pronounced upon cooling. The position of the second peak changes monotonously: it moves to higher Q values on decreasing temperature over the entire composition regime.
its position moves slightly to lower Q values with decreasing temperatures in mixtures with low methanol content, then stays nearly stable around the equimolar mixture, and later moves to higher Q values in the methanol rich compositions.
The X-ray and neutron weighted total scattering structure factors (F X (Q) and F N (Q), respectively),

Partial radial distribution functions
As simulated total structure factors reproduce measured data satisfactorily, and thus, the simulations are validated against available diffraction data, PRDFs obtained from MD calculations may be considered as representative. Note that since the number of atom types in the system is 6 (methanol: C, O, methyl H denoted as H, hydroxyl H denoted as HM; water: OW and HW), the number of partial radial distribution functions is 21. That is, it would be impossible to separate the partials, using only the two experimental datasets, via any other methods.
Temperature dependence of the PRDFs related to H-bonding is shown for xM = 0.7 in Figure  2xM+4(1-xM) (this is best observable for methanol concentrations xM ≥ 0.5).
The number of H-bonds between different kinds of pairs of molecules is presented in Figure 8. Based on the above findings, we conjecture that lowering the temperature forces enhanced mixing of the two species in certain concentration and, perhaps more importantly, temperature regions. It is important to point out that below 193 K, the number of like-like (particularly of water-water) HB-s start to rise again on further cooling. It may be speculated that this connected with the precursor of water freezing out of the liquid mixtures in these regions. This may be an explanation for that in some of the mixtures, the number of HB-s between unlike molecules start to decrease below ca. 193 K (cf. Fig. 8 parts b) and d) ). Note, however, that statistical uncertainties, as well as simulated annealing related issues in MD simulations, would necessitate at least one order of magnitude longer simulation times for confirming these speculations.
According to the tendencies mentioned above, taking by and large, three concentration regions may be identified, with respect to changes of hydrogen bonding on cooling. At low methanol concentrations, in the ca. xM ≤ 0.3 region, the number of water-water and water-methanol pairs are increasing similarly with decreasing temperature, while the number of methanol-methanol pairs is nearly constant. (Similar trends were observed in Ref. [8].) In the 0.3 < xM ≤ 0.7 region the number of waterwater pairs is growing slower than that of water-methanol pairs, whereas the number of methanolmethanol pairs is constant (or slightly decreasing). In the methanol rich region, xM > 0.7, the number of water-water pairs is nearly constant, the number of water-methanol pairs is growing as in the other regions, and the number of methanol-methanol pairs is also increasing with decreasing temperature.
The number of H-bonded molecules, water and methanol together, around both water and methanol is larger at lower temperatures (NW and NM, Figs. 8c and 8f). However, the behavior of the two components at the lowest investigated temperatures is different: the number of H-bonded molecules around methanol is saturating for mixtures where xM > 0.4, but the number of H-bonded molecules around water is increasing significantly even at the lowest investigated temperatures in each mixture. This may be explained by remembering that water molecules are easily capable of forming four hydrogen bonds, whereas methanol molecules, due to simple steric considerations, are hardly capable for more than two in practice.
The concentration dependence of the number of H-bonded pairs at some selected temperatures are shown in Figure 9. The number of water molecules around both water and methanol (NWW and NMW, Figs. 9a and 9c) is increasing as methanol concentration decreases. The shape of the curves is not linear: they both differ from linear upward, i.e. the number of hydrogen bonds to water molecules is higher than it would be proportional to the concentration. At 300 K the number of H-bonded methanol molecules around water and methanol (NWM and NMM, Figs. 9b and 9d) also show non-linear concentration dependence, but the deviation from linearity is negative for these curves. The same behavior was found at 300 K previously [14] by MD simulations using the SPC/E water model. As temperature is decreasing, the shape of the curves changes slightly, as the number of H-bonds is increasing. The most significant change can be observed in terms of the number of methanol molecules around water (NWM, Fig. 9b): as temperature is decreasing the deviation from linearity becomes at first less negative, and at 193 K it even switches to positive. These observations suggest that (at least at high methanol concentrations) preference of the H-bonded mixed pairs is increasing with decreasing temperature. This is in line with the observation made above, i.e., that decreasing temperature enhances mixing in certain concentration ranges.

Conclusions
Synchrotron X-ray and neutron diffraction experiments have been conducted, as a function of decreasing temperature, over the entire composition range of methanol-water liquid mixtures. Molecular dynamics simulations have been made use of for interpreting experimental data. The all atom OPLS/AA force field model for methanol has been combined with both the SPC/E and TIP4P/2005 water potentials. Based on results obtained from these investigations, the following specific statements can be made: The TIP4P/2005 water model was found to be somewhat more successful in reproducing measured X-ray and neutron diffraction data in the reciprocal space.
(ii) Hydrogen bond related partial radial distribution functions, as calculated from simulation trajectories, show sharpening first (and, wherever well distinguishable, second) maxima and minima. The positions of extrema do not change significantly.
(iii) As a general trend, the average number of hydrogen bonds increases upon cooling. This is strictly true for the 'water-all' and 'methanol-all' cases.
(iv) Interestingly, the number of hydrogen bonds between methanol molecules slightly decreases with lowering temperatures in the concentration range between ca. 30 and 60 mol % alcohol content. The same is valid for water-water hydrogen bonds above 70 mol % of methanol content, from room temperature down to 193 K.
(v) The effects of decreasing temperature is most significant for hydrogen bonding between unlike (water and methanol) molecules; this indicates that decreasing temperature enhances mixing between the constituents.
[        The all atom OPLS-AA [3] force field was applied for methanol, whilst for water two models, the SPC/E [4] and the TIP4P/2005 [5] were tested. Non-bonded interactions were described by the 12-6 Lennard-Jones interaction and the Coulomb potential (see Eq. 1): where rij is the distance between particles i and j, qi and qj are the partial charges on these particles, ε0 is the vacuum permittivity, and εij and σij represent the energy and distance parameters of the LJ potential. The LJ parameters (εii and σii) and the partial charges applied (qi) to the different atoms are collected in Table S1. The εij and σij parameters between unlike atoms are calculated as the geometric average of the homoatomic parameters (geometric combination rule, in accordance with the OPLS/AA force field).
Intramolecular non-bonded interactions between first and second neighbors were neglected, whereas between third neighbors (atoms separated by 3 bonds, HC -HO bonds) they were reduced by a factor of 2. The intramolecular (or bonded) forces considered here are the bond-stretching (2body), angle bending (3-body) and the dihedral angle torsion (4-body) interactions. Bond lengths in methanol molecules were fixed using the LINCS [6] algorithm, while bond angles and torsional angles were flexible. The rigid water geometry was handled by the SETTLE algorithm [7]. Bond lengths, equilibrium angles and force constants are given in Table S2.
The smoothed particle-mesh Ewald (SPME) method [8,9] was used for treating Coulomb interactions, using a 20 Å cutoff (15 Å for methanol concentrations 0.1 and 0.2 molar fraction, due to the smaller box sizes) in real space. Non-bonded LJ interactions were cut-off at 20 Å (15 Å for methanol concentration of 10 and 20 mol%), with added long-range corrections to energy and pressure [10].
Initial configurations for the NPT simulations at T = 300 K were obtained by placing the molecules into the simulation box randomly, following an energy minimization using the steepestdescent method (random configuration method). At lower temperatures the final configuration of the previous temperature point was used as initial configuration. The equations of motion were integrated via the leapfrog algorithm, the time step was 2 fs. At each temperature, at first a short (0.2 ns) NVT run was performed using the Berendsen thermostat [11], with τT = 0.1, for relaxing the system to the target temperature. After that, the NPT ensemble was used. The temperature was kept constant by the Nose-Hoover thermostat [12,13], with τT = 2.0, while the pressure was kept at p = 10 5 Pa, by the Parrinello-Rahman barostat [14,15], using a coupling constant of τp = 2.0. After a 2 ns equilibration period, a 2 ns production run was completed, from which the densities were calculated.
Initial configurations for the NVT simulations were either obtained by the random configuration method, or were adopted from the corresponding NPT simulations. (The latter method was used for low temperatures and high methanol content mixtures, in order to avoid artifacts resulting from close packing and low mobility of molecules). The leap-frog algorithm was used again, with the same time step as for NPT runs (2 fs). Two equilibration steps were performed before the production run: during the first, short one (0.2 ns), the Berendsen thermostat was used; after that, and also for the production run, the Nose-Hoover thermostat was activated, with the same coupling constants as before.
Trajectories were saved in every 200 ps, for the duration of 20 ns: this way, 101 configurations could be used for further analyzes. Partial radial distribution functions (PRDF, gij(r)) were calculated from the collected configurations, by the 'gmx_rdf' programme of the GROMACS software. The model structure factor can be obtained from the PRDFs, according to the Faber-Ziman formalism [16], by the following equations: where Q is the amplitude of the scattering vector, and ρ0 is the average number density.
The XRD (F X (Q)) and ND (F N (Q)) total structure factors can be composed from the partial structure factors Sij(Q) as: and where wij X,N denotes the X-ray and neutron scattering weights. For X-rays it is given by equation (5): Here δij is the Kronecker delta, ci denotes atomic concentrations, fi(Q) is the atomic form factor.
The neutron weight factors are: where bi is the coherent neutron scattering length.
The total structure factors obtained from simulations were compared with the measured curves by calculating the goodness-of-fit (R-factor) values: where Qi denote the experimental points, 'mod' indicates the simulated and 'exp' the experimental curves.