The observation of oxygen—oxygen interactions in ice

In this article, we present a series of first principle density functional theory simulation results which have, reproduced the vibrational spectra for ice Ih and VIII obtained by neutron scattering. Detailed analysis of the pressure-dependent simulation results shows that there is an oxygen–oxygen (O–O) interaction between the nearest water molecules participating in the H-bonding and this interaction is 1.5 times stronger than the H-bond at an O–O distance Roo of 2.7 Å observed in ice I. The O–O interaction strength increases rapidly to 2.3 times the H-bond at Roo = 2.5 Å. Analysis of the simulation results for ice VIII, shows that the strong O–O interaction only exists between H-bonded water molecules (i.e. O-H—O) and not for non-H-bonded water molecules, even when the O—O distances are similar. This discovery completely changed our view of H-bonding in water and ice, redefining the term of hydrogen-bonding between water–water molecules. Incorporation of this term may also help us to explain water anomalies. In addition short H-bonds (Roo < 2.7 Å) exist in DNA and proteins, and are thought to play a role crucial in enzyme catalysis and function. Hence this result could have major implications in our understanding of enzyme-catalysis and biology.


Introduction
An accurate description of the H-bonding between water-water molecules is widely recognised as the crucial factor in the understanding of the properties of water (water anomalies in particular). Despite a large number of water potentials available in the literature, from simple three points charge model such as MCY [1] and SPC [2], to the four points model TIP4P [3] and more recent five points TIP5P [4], none are able to reproduce all water anomalies and in particular the vibrational phonon density of states (PDOS) obtained from inelastic neutron scattering (INS) for the exotic ices [5,6]. Hence there is still neither consensus nor model which can be widely applied. One important reason for failure is the lack of the flexibility of the covalent bond O-H which could change to respond the hydrogen bonding (O-H) changes. Recently Sun et al suggested that such changes could be significant in altering the properties of water and ice [7]. In fact, there were a number of attempts in literature to represent the flexibility of the covalent bonding in water/ice, such as the flexible-SCP [8] and TIP4PF [9] by adding a spring force on the covalent bond, this simple harmonic repetition of the intra-molecules motion is of limit success [10]. Kumagai et al has proposed a flexible potential (named KKY potential) which has three separate terms: V HO (r), V HH (r) and V OO (r) in order to account the intra-and inter-molecular interactions separately [11], our MD simulation for the INS spectra based on the KKY potential was also unsuccessful [12].
In order to understand the relationship between intra-and inter-molecular motions and to accurately describe the coupling the covalent and H-bonding in water, ab initio quantum mechanical simulation has demonstrated the best way forward [13][14][15]. In our previous studies [16,17], we have applied first-principle density function theory (DFT) to reproduce the INS spectra obtained and the simulation results show a good agreement with the measured data, implying the ab initio technique used is very credible in producing the vibrational PDOS for ice Ih without need of the empirical input. This was due to the newly updated CASTEP [18]  has made structural and energy convergences much precise and effective for H-bonded systems and hence has enabled us to reproduce the PDOS measured by INS. In this article, we present further analysis of these simulation results with intention to understand the force field generated by the ab initio DFT for the subsequent PDOS calculations in order to compare with the force constants used in the early model [19].

Simulation methods
Using ab initio DFT, a series of the simulations for different formats of exotic ice crystal structures were performed. The generalised gradient approximation (GGA) function of PW91, HCTH and RPBE were used to describe the exchange-correlation effects (results show little differences among the different functions and hence the PW91 results were presented here). The geometry optimization calculations for the ice structures were first conducted. In order to get refined results we set higher restrictions than the fine-quality by default, the convergence tolerance of the energy is set up to 1.0E-6 eV, the maximum force is 1.0E-5 eV, and maximum displacement is 1.0E-4 Å. The maximum iteration is expanded to 1000 steps to match the accuracy above, the maximum step size is 0.1 Å. In addition, a bigger plane wave energy cutoff of 800 eV for the ultrasoft pseuedopotential has been used with better results. The van der Waals force was also examed by adopting the DFT-D in the CASTEP code, its effect to the stability of the ice structures has reported [20], a small effect of redshifts (∼2 meV) for the main peaks of the vibrational spectra was observed. The PDOS were calculated using the CASTEP module with finite displacement method [21]. The results were compared with experimental PDOS obtained by INS [5]. Further analyzing the simulation results, the force constants produced by the CASTEP for the PDOS calculations were obtained from the output files.

The result and discussion
The success of the simulation method can now reproduce PDOS not only for ice Ih, but also for other exotic phases of ice, such as ice VIII as shown in figure 1, including the two optical peaks at 28 and 37 meV in ice Ih. This result suggests that the 2 optical peaks can be generated from the simulation without invoking the two H-bond forces proposed in the past [19]. Slight shift of the two optical peaks was found in the simulated spectrum of ice Ih towards higher energies, at 32 and 41 meV, indicating the PW91 functional has overestimated the strength of the H-bonding. Other functionals tested also show the blue shifts with different degrees without change of main features of the PDOS spectra. Adding DFT-D in the simulations shows a small improvement of the blue shifts without change of the main features of the spectra. The librational band of ice Ih was also blue shifted about 10 meV. For ice VIII, the simulation reproduces the main features of INS spectrum (with similar degree of blue shifts), such as the three peaks in the translational region (<50 meV) and the two sharp peaks in the left hand side of the librational region at about 65 meV. The small peak in the right hand side of librational region at about 100 meV was not shown in the INS spectra due to large Debye-Waller effect in the measured spectra which smeared the spectra dramatically at higher energy transfer. Using different INS spectrometer at ISIS, we can see this peak clearly [5].
The interesting results obtained from these simulations are for the pressurised ice Ic structure, as shown in figure 2(a) at a range of different water-water separation distances (R OO from 2.1 to 3.5 Å). We chose the cubic structure of ice Ic instead of ice Ih because ice Ic has a higher symmetry and fewer water molecules in the unit cell than ice Ih, but its PDOS and INS spectrum are identical to ice Ih (and same as the spectrum of low density of amorphous ice [6]). For each Roo distance, the corresponding lattice constants were fixed during the simulation and the O and H atoms allowed to move, aiming to minimise the total energy of the structure.
After the structural minimisation, the CASTEP calculates PDOS by applying a small displacement of each atom around their equilibrium positions in the ice structure along the x-, yand z-directions, small energy variations along these directions were obtained. The force for the atom i, F i (=dV/dr i ), is simply the derivative of the total energy V. By applying a second small displacement for the atom j, the force constant K ij (=d 2 V/dr i dr j ) is obtained for the pair atoms i and j. A force matrix [K] for the unit cell which is a 3Nx3N matrix (where N is total number of atoms in the unit cell) were constructed based on above procedure, hence PDOS can be calculated by a subroutine in the CASTEP program. The structure of water at the O-O distance at 3.5 Å (corresponding to a negative external pressure of −1.2 GPa) is almost like a gas phase, i.e. the H-bonding at such a large separation distance is negligible. The spectrum for this structure can be considered as a group of almost isolated water molecules with very weak intermolecular interactions. As one can see from the top curve in figure 3, the translational and rotational modes shifted to lower energies (with translational modes cut-off at 31 meV, in comparison of 42 meV for the normal ice) reflect the lack of strength of the H-bonding. When R OO decreases, the translational and librational frequencies increase. The intramolecular modes at frequencies above 400 meV are two delocalised modes of the symmetric (ν 1 ) and anti-symmetric (ν 3 ) O-H stretching which are well separated from each other as seen experimentally for isolated water molecules. When the O-O distance is further reduced, the spectra gradually develop into that of ice Ih as shown in figure 1. At the R OO distance of 2.7 Å the spectrum is more or less identical to that of ice Ih at ambient pressure as measured by INS. With O-O distances less than this (equivalent to the structure under very high external pressures), the translational modes increase towards higher energies, as do the librational modes. For the intramolecular vibration region, the ν 1 and ν 3 modes decrease, with the whole stretching band spreading further. For R OO at 2.3 Å (corresponding to an external pressure of 27 GPa) the H-bond length equals that of the internal covalent bond at 1.15 Å, i.e., the H atom actually sits in the middle of the two nearest oxygen atoms as shown in figure 2(b). Technically, the molecular character of the ice crystal breaks down at this R OO and its structure becomes symmetric, as observed experimentally [22].
After detailed analysis of these force field from the output of the ab initio simulations which is an intermediate stage of the PDOS calculation, all pairs of force constants in the force matrix were obtained. The significant ones, such as the covalent bond force constant K CB , H-bond force constant K HB were plotted in figure 4. As one can see, at large R OO of 3.5 Å, K HB is small (about 0.53 eV Å −2 ) indicating a weak H-bonding, and K CB (=47.7 eV Å −2 ) is very large. When the O-O distance decreases the K HB increases and the K CB decreases, at R OO =2.7 Å (i.e. at ambient pressure) the K HB and K CB are 2.21 and 34.1 eV Å −2 , respectively, and they are almost identical to the values used in our earlier model [19], indicating the data obtained in the ab initio simulations are credible. When the structure becomes symmetrical at R OO =2.3 Å the force constants K CB and K HB are equal at 5.4 eV Å −2 as expected. The error bars for both the K CB and K HB (obtained from the spread of the force constants on the different bonds in the structure or divergent from the averaged values) are too small (less 1%) to plot in the figure 4, indicating that the structure is well relaxed. Furthermore, no second H-bond force constant presents in the force field matrix indicating only one H-bonding strength in the ice structure.
The most interesting result of the simulation was a discovery of an additional O-O interaction, K OO , between adjacent oxygens which becomes pronounced at R OO less than 3.1 Å with the K OO :K HB ratio of 1:1, but increases rapidly when the R OO decreases (see table 1). In fact this force constant increases to 2.3 times of the K HB at R OO =2.5 Å, indicating a very strong interaction between the two O atoms. This additional O-O interaction on top of the H-bonding has not been considered in the past however it plays an important role in reproducing the INS spectra and it has a completely different R-dependence. This additional contribution may shine a new light on our understanding of H-bonding in water, and could be possible to explain why water properties become abnormal in liquid or other phases when the R OO is less than 3 Å.
Extending the analysis further, we also applied the same considerations in pressure-dependent simulations for ice VIII [17] which has a cubic structure containing two sets of interpenetrated ice Ic sub-lattices as shown in figure 2(c) and hence each water molecule (or oxygen atom) has 8 nearest neighbours, 4 of them are H-bonded

Conclusion
In this article, we illustrated the effectiveness of current ab initio techniques when it is dealing with H-bonded systems. The technique has successfully reproduced the PDOS for ice Ih, Ic, VIII and other structures we simulated, in particular the two optical modes at 28 and 37 meV which have previously proven recalcitrant to simulation. Detailed analysis of the force constants obtained from the ab initio simulations shows there is only one H-bonding strength found in the structure with a very small spreading (less than 1%). As we indicated in our original model [19], two forces are an essential element for the reproduction of the two peaks found in the PDOS obtained by INS. However this requirement is met by the observation of the additional K OO interaction between H-bonded O atoms. More importantly, this additional force could provide the necessary mechanism to explain a range of water anomalies, for instance the complex phase diagram, i.e. its morphism, this is probably because the HBs become very easy to buckle under the extra strong O-O interaction if the O-H-O is not in straight line. This O-O interaction, K OO has not been realised and observed experimentally in the past, perhaps because it is easily concealed by the K HB and K CB . As we know, the force constants are the second derivative of the interaction potential, i.e. K ij =d 2 V/dr i dr j . Most classical rigid water-water potentials V(r) have a L-J term and charge term (possibly with additional polarization term, pol) representing the H-bond: where the distance of r is the distance of O-O atoms, i.e. r=R OO and the r' is the distance of charged Oand H + . Based on the function, we will be able to obtain the force constants of K HB to represent the H-bonding. The intra-molecular interaction can adopt a spring constant as provided by the flexible-SPC or TIP4PF. However our result shows that apart from the K HB there is additional K OO which cannot be obtained based on equation (1) at the same time, in other words, we cannot obtain two separate values of force constant at a fixed r with the same local geometry (i.e. the charge arrangement, etc) and hence the curves for K HB and K OO shown in figures 4 and 5 cannot be obtained based on the current forms of classic potentials [1][2][3][4] without introducing an additional term in their formula. In addition we found that the strong O-O interaction only occurs between the H-bonded O-O atoms, indicating that the H atoms between the two O atoms plays an essential role in bringing the two O together, a detailed mechanism is still unclear. Considering that short H-bonds <2.7 Å are common in DNA and proteins [23], and have been observed clustered in the catalytic sites of enzymes [24], this result suggest that the strong O-O interaction could also play a very important role in the understanding of the efficiency and specificity of biological catalysis and the evolution of enzymes.