Decreasing temperature enhances the formation of sixfold hydrogen bonded rings in water-rich water-methanol mixtures

The evolution of the structure of liquid water-methanol mixtures as a function of temperature has been studied by molecular dynamics simulations, with a focus on hydrogen bonding. The combination of the OPLS-AA (all atom) potential model of methanol and the widely used SPC/E water model has provided excellent agreement with measured X-ray diffraction data over the temperature range between 298 and 213 K, for mixtures with methanol molar fractions of 0.2, 0.3 and 0.4. Hydrogen bonds (HB-s) have been identified via a combined geometric/energetic, as well as via a purely geometric definition. The number of recognizable hydrogen bonded ring structures in some cases doubles while lowering the temperature from 298 to 213 K; the number of sixfold rings increases most significantly. An evolution towards the structure of hexagonal ice, that contains only sixfold hydrogen bonded rings, has thus been detected on cooling water-methanol mixtures.

Despite the voluminous literature available on liquid methanol and methanol-water mixtures, the flow of new investigations on the subject is continous (see, e.g., refs 1-5). The reason, apart from the enormous industrial significance of methanol, is that the molecules of methyl-alcohol, CH 3 OH, are the closest organic analogues, by their size and shape, of water molecules, but they contain a hydrophobic methyl (CH 3 -) group. These ambivalent properties of methanol molecules give rise to a complex structure of its mixtures with water: hydrogen bonds (HB-s) may form between water-water and methanol-methanol, as well as between water-methanol and methanol-water pairs. Characteristics of hydrogen bonds in methanol and water-methanol mixtures have been considered by numerous publications [5][6][7][8] ; for mixtures, these analyses have often been conducted as a function of composition 1, 2, 9 . Interestingly, hydrogen bonding as a function of temperature has not yet been analysed in detail in water-methanol mixtures. One of the reasons may be the lack of suitable experimental data: while diffraction results for ambient conditions are abundant, both by X-rays 2, 10-12 and by neutrons 6,7 , to our best knowledge, only one X-ray data set is available for temperatures down to the freezing point 12 . This lack of information has turned our attention to structural investigations at lower temperatures.
In a recent work, some of the present authors reported a detailed analysis of molecular dynamics computer simulation data for the microscopic structure of water-methanol mixtures in the entire composition range and performed comparisons with the experimental results obtained by X-ray diffraction 2 . The principal focus in that publication was on the changes of the experimental and theoretical total structure factors of mixtures as a function of composition, considered at room temperature and ambient pressure. On the molecular dynamics (MD) computer simulation side, we analyzed the SPC/E 13 and TIP4P-Ew 14 water models, in combination with the OPLS-AA methanol model 15 .
Encouraged by an overall satisfactory performance of the molecular dynamics simulation data for describing experimental trends, here we take temperature dependent X-ray diffraction data of Takamuku et al. 12 and perform molecular dynamics simulations for compositions and temperatures of the diffraction experiments, using the same interatomic potential functions (SPC/E for water and OPLS-AA for methanol) as previously 2 . The resulting particle configurations have been analyzed by the tools reported in refs 4 and 9 for identifying hydrogen bonds and find topological characteristics of their networks.
Ring structures composed by hydrogen bonded molecules are an easily recognizable type of hydrogen bonded assemblies. Their significance is indicated by the fact that the structure of one of best known material in everyday life, crystalline ice, is entirely built of rings that contain 6 hydrogen bonded water molecules 16 . Cyclic entities have previously been found in simulated water-methanol mixtures at room temperature over the entire composition range 9 , although the number of rings decreased rapidly as methanol concentration grew. In this work, the temperature evolution of the number of cyclic entities, that are kept together via hydrogen bonds, is followed. We expect that the insights provided by the present investigation will contribute significantly to our understanding of structural changes that occur in hydrogen bonded systems upon lowering the temperature.

Results
Total structure factors. Figure 1 compares X-ray weighted total structure factors, S(Q), as measured by X-ray diffraction 12 and as calculated from the present molecular dynamics computer simulations for the mixture with Xm = 0.4 (Xm is the molar fraction of methanol). The applied interatomic potentials have proven to be adequate for reproducing X-ray diffraction data over the temperature range between 298 and 213 K, and for methanol molar fractions of 0.2, 0.3 (not shown) and 0.4. Simulated curves follow the same trends with temperature as those obtained from diffraction experiments; moreover, calculated S(Q)-s match the corresponding measured functions at an at least semi-quantitative level. Deviations between experiment and simulation could only be observed around the second maximum of the total scattering structure factor, similarly to what had been found for room temperature diffraction data 2 . The very good agreement between diffraction data and simulated total scattering structure factors provides a strong basis for detailed structural analyses. In other words, the present calculations may be considered as validated from the structural point of view, on the basis of the agreement with X-ray diffraction experiments.
Radial distribution functions. Partial radial distribution functions, PRDF, are inaccessible from diffraction experiments for water-methanol mixtures and therefore, only computer simulation results may be presented in Fig. 2. Possessing atomic coordinates, it would be possible to calculate all the PRDF-s; here we restrict ourselves to displaying the functions that are related to hydrogen bonding, i.e., all possible combinations of O-O and O-H partials. Each of the 7 PRDF-s show the same trends as temperature decreases: maxima and minima become slightly, but systematically sharper. Also, slight shifts in terms of the positions of first O-O maxima and minima can be observed, due to the slight change of the density with decreasing temperature. What is important to put on the record at this stage is that the observed changes in terms of two-body correlations (the PRDF-s) are consistent with common expectations and that these changes are small. These findings are in qualitative agreement with the suggestions of Takamuku et al. 12 .
Hydrogen bond characteristics. Hydrogen bonds have been identified both by energetic and/or geometric considerations, as in previous publications 2,4,9 . According to these preliminaries, two molecules are identified as H-bonded if the non-bonding distance between an oxygen and a hydrogen atom is less than 2.5 Å, and (a) the O…O-H angle is smaller than 30 degrees (geometric definition), or (b) the interaction energy between the H-bonded molecules is less than −3.0 kcal/mol (energetic definition). All calculations concerning H-bonding and the H-bonded network have been performed by the HBTOPOLOGY code, originally developed by Chihaia et al. 17 .
Based on Table 1 the following observations can be made: (i) The number of energetically defined hydrogen bonds is lower, by about 5 to 10%, than the number of hydrogen bonds identified via purely geometrical means. This difference is roughly constant over the entire composition and temperature range considered here, so that all the trends that may be identified for one definition is valid for the other one, as well. Nevertheless, key quantities will be mentioned for both definitions, so that it is ensured that the main findings are not influenced by the H-bond definition. The (huge) extent of the increase in terms of the number of cyclic entities would be impossible to predict on the sole basis of points (i) to (iv); it is therefore thought that this finding may termed as 'surprising' and 'unexpected' . In contrast, data concerning the number of H-bonds in Table 1  The last finding of shown in the above list ['(v)'] led us to look at primitive cyclic entities (the ones that cannot be split into smaller rings 17 ) more closely. Figure 3 displays ring size distributions, as a function of temperature, for the mixtures simulated in the present study, for H-bonded cycles containing up to 8 molecules.
The first thing to notice is that for both definitions of H-bonds and for each composition the overall number of cycles becomes progressively larger as the temperature approaches the experimental freezing temperature. The increase is significant for 5-, 6-and 7-membered rings, and it is outstanding for the 6-membered cycles. At higher temperatures and higher methanol contents, 5 membered rings are slightly more abundant than 6-membered ones, but at the lowest temperature, 6-membered primitive cycles become dominant at each concentration. The composition of 6-membered rings is, apparently, not sensitive to decreasing the temperature: for instance, for the X m = 0.2 mixture at 298, as well as at 243 K, 58% of such cyclic entities contain no methanol molecule at all, whereas 1 and 2 methanol molecules can be found in 36 and 5%, respectively, of the 6-membered rings. The compositions of 5-and 7-membered cycles are similarly stable against temperature.
The fact that 6-membered rings contain significantly more water molecules than it would follow from the composition of the mixtures indicates that these homogeneous systems look microscopically somewhat heterogeneous, particularly as temperature goes down -it might even be said that one component (6-membered rings of water molecules) is dispersed in a homogeneous 'matrix' . This point would certainly be important to investigate for mixtures with higher alcohol contents.
In order to give an idea about the populations of 6-membered rings in the X m = 0.2 mixture at 298 and at 243 K, Fig. 4 draws these cycles at these two temperatures (these drawings are taken from typical single particle configurations, using the energetic definition of H-bonds; note that for the geometric definition the number of sixfold rings is even larger). While the presence of 6-membered rings at room temperature is already significant, they seem to fill the entire simulation box at the lowest temperature considered (243 K). The change is dramatic; we make an attempt of providing a simple, qualitative explanation in the next section.

Discussion
In the light of the above findings, it is rather unfortunate that experimental data are available for only these compositions, all in the water-rich regime. This is the reason why we are unable to extend experimental data based investigations to methanol-rich compositions. For pure water (0 molar % methanol) it is known 18, 19 that the SPC/E potential model cannot reproduce diffraction data at room temperature; for this reason, pure water had to be left out from this study.
To the best of our knowledge, such greatly enhanced ring formation in water-alcohol mixtures on cooling has not been spotted in simulated structures that were in good (i.e., at least, semi-quantitative) agreement with measured diffraction data. For instance, the original paper of Takamuku et al. 12 suggests that 'alcohol-water mixtures are frozen with keeping the structure of the dominant clusters formed in the mixtures at 25 °C' , which clearly is different from what has been found here. The extent of the increase in terms of the number of H-bonded rings would be impossible to predict on the sole basis of the increase of H-bonds. Based on the present (although admittedly, limited) sample it may be conjectured that this is an essential phenomenon that may well be worth monitoring in other hydrogen-bonded systems, as well.
We suggest that a handwaving, qualitative explanation for the striking increase of cyclic entities could be that in hydrogen-bonded systems with substantial amounts of water, joining 'dangling' hydrogen bonds is an effective way of decreasing the entropy of the system. The presence of water is thought to be important since in all the known crystalline forms of water, H-bonded rings dominate; for instance, the most common form of crystalline water, ice Ih, consists of 6-membered rings exclusively. As pointed out in the previous section, primitive   cycles identified in our simulated structures contain mostly water, with only about one methanol molecule per 6-membered ring on average. The clear domination of 6-membered rings in each low temperature particle configuration suggests that the appearance of a large number of cyclic entities is, while certainly lowering the configurational entropy of the system, likely to be a precursor of freezing. As a final thought, it may be concluded that even though the enhanced ring formation might have looked extraordinary at first sight, it is quite easy to find ordinary, simple reasons for why it should happen in water-dominated alcohol-water mixtures. This is why it is interesting why this phenomenon has not been spotted earlier; one reason is a lack of suitable diffraction results to compare with -and it would be a desirable line of investigations to extend our knowledge base in this respect. Studying pure water in a similar fashion would be imperative, in order to check what happens in the pure substance that is known to form 6-membered rings in its common crystalline phase ice Ih. In parallel, it would be intriguing to explore the structure of pure alcohols as a function of temperature, since their crystalline phases do not consist of rings of hydrogen bonded molecules.

Methods
All the simulations of water-methanol model mixtures in this work were performed by using the GROMACS molecular dynamics package 20 (version 5.1.1) in the canonical ensemble at each temperature (298 down to 213 K) and at ambient pressure.
We chose the SPC/E 13 model for water, combined with the OPLS-AA model 15 for methanol, as this combination was shown to work well for methanol-water mixtures at room temperature 2 . The geometric combination rules, ε ε ε = ( ) ij ii jj 1/2 and σ σσ = ( ) ij ii jj 1/2 for the interaction between unlike species (cross, non-bonded interactions) were used, together with the same rules for pure components. Non-bonded interactions were cut-off at 11 Å. In all the cases, bonds within methanol molecules were maintained fixed, whereas the angles and torsional angles were flexible. For the sake of maintaining the geometry of the water molecules the LINCS algorithm 21 was used. The Particle Mesh Ewald method was used for treating Coulomb interactions 22 ; the cutoff was again 11 Å.
Initially, particles were placed in a cubic simulation box randomly. Periodic boundary conditions were used. The number of molecules of both species in the simulation box for the three compositions is given in Table 2. Berendsen thermostat with τ T = 0.1 ps was applied 23  After equilibration, production runs were conducted. Simulations last 10 ns after equilibration. For hydrogen bond analyses, samples of 101 particle configurations, 10 ps apart, were collected.
For the identification of cyclic entities, the computer programme 'HBTOPOLOGY' 17 was applied; further applications of this code can be found in refs 4 and 9. We note that quite a few publications report on algorithms that can be used for the same purpose (see, e.g., refs 24-27 and references therein); HBTOPOLOGY was chosen since the source code was kindly made available to one of the present authors (I.B.) and therefore the transferability between various hydrogen bonded systems could be ensured.
The algorithm starts with setting up the connectivity list, according to the actual definition of hydrogen bonds. Each oxygen atom is taken as a 'node' in the network. The ring search is carried out on the basis of the connectivity of the nodes, considering one after another each node as starting-node. For a cyclic entity, the starting node has to be identical to the end node. The programme simply explores all possible paths that return to the starting node while touching a limited maximum number of nodes along the path. This pre-defined number was 6 in the original publication 17 ; in the present case, cycles of up to 8 nodes have been searched for. Results presented in Fig. 3 concern primitive rings, i.e. ones that cannot be decomposed into smaller rings.   Table 2. Number of molecules of water and methanol in MD simulations, molar fraction of methanol and weight concentration.