Hydrogen-bond potential for ice VIII-X phase transition

Repulsive force between the O-H bonding electrons and the O:H nonbonding pair within hydrogen bond (O-H:O) is an often overlooked interaction which dictates the extraordinary recoverability and sensitivity of water and ice. Here, we present a potential model for this hidden force opposing ice compression of ice VIII-X phase transition based on the density functional theory (DFT) and neutron scattering observations. We consider the H-O bond covalent force, the O:H nonbond dispersion force, and the hidden force to approach equilibrium under compression. Due to the charge polarization within the O:H-O bond, the curvatures of the H-O bond and the O:H nonbond potentials show opposite sign before transition, resulting in the asymmetric relaxation of H-O and O:H (O:H contraction and H-O elongation) and the H+ proton centralization towards phase X. When cross the VIII-X phase boundary, both H-O and O:H contract slightly. The potential model reproduces the VIII-X phase transition as observed in experiment. Development of the potential model may provide a choice for further calculations of water anomalies.

the decrease of the molecular coordination, H-O covalence bond contracts spontaneously to lower the cohesive energy, and hence O:H nonbond is polarized, leading to the surface lubricity 35,37 . This force is beyond the conventional intra-and inter-molecular interactions but depending on the existence of the O:H-O link, i.e. if there is no O:H-O bond like between two O atoms in difference sublattice in ice VIII, hidden force (3 times stronger than Van de Vaals force) interaction would not show 38 . However, it has still seldom reported to address the hidden force of O:H-O bond in a H-bond model, although the force plays a significant role in altering the properties of water and ice.
In this paper, in order to investigate the strange behaviour of proton centring in phase VIII-X transition, we clarified a hidden force model considering the quantum interactions between electron clouds of covalent bond and nonbond and reproduced the anomalous behaviour of proton in ice VIII-X transition.

Potential Model Proposed
Based on our extensive investigation of water anomalies 21,34-36 , we proposed hidden force model for dynamics of the "O: H-O" bond under external pressures, as shown in Fig. 1a. In the model, "O: H-O" bond is similar to be connected by three "non-harmonic converted to anharmonic springs". Black spring represents covalent bonding of H-O (denoted as subscript C); Grey spring represents van der Waals interaction between lone pair (blue circles) and proton of O:H (denoted as subscript H); Blue spring represents repulsive interaction between bonding electron pair and lone pair (denoted as subscript OO, but the interaction is not the conventional O-O coulomb repulsion). A water molecule is surrounded by four nearest neighbours (Pauling's ice rule). Four identical O: H-O forms a tetrahedron as the basic unit of water structure ( Supplementary Fig. S1 where proton disordered symmetrization occurs with delocalized proton in a relatively broad potential well. The delocalization of proton (or proton-transfer) may attributed to quantum effect of nuclei 16,39 , or to thermal fluctuations obtained by AIMD 15 . In this work, we focus on the proton centring due to electron interactions within the O:H-O bond.
Without external force, taking centre proton as reference, three springs can reach to equilibrium with their forces along the directions as indicated in Fig. 1a. Due to the f oo , both O atoms are pushed away a little from their lowest-energy positions in the potentials of isolated bonds. Thus, the initial potential recovery forces f C and f H will both point inward to center proton; while repulsive force f OO will point outward. Supposing ≅ + x x x OO H C , at equilibrium, total forces (force values) added on O atoms should both be zero: where, x C0 and x H0 (length values in positive numbers) are the equilibrium position of O atoms counting form x = 0. Thus, in equilibrium without external force, Considering an identical force f p at both sides, the forces satisfy: , combining eq. (2) and eq. (3), we can get: Eq. (4) tells us: i) The values of δx C and δx H should be determined by the potential forms of V H and V C . Since O:H bond is much longer and weaker than H-O bond, V H is more flat (smaller curvature) than V C , i.e. |k C | > | k H |. Thus, the value of δx H is always larger than δx C . ii) If δx C and δx H have opposite signs, i.e. x C and x H do not both increase or decrease together, k H and k C must have opposite sign with each other, which means the curvatures of potentials are opposite, as shown in Fig. 1a. iii) If δx C and δx H have same signs, i.e. x C and x H both increase or decrease together, k H and k C must have the same sign.
ii) Stage II (ice X): external force is strong enough to make x H shrink over the curvature turning point (star in Fig. 1) and k C • k H > 0. External force contract both segments, pushing O atoms back towards the energy-lowest positions, compressing all the three springs. In this case, f P and energy increases for both segments provides the sharp increase of the V OO :

Calculations and Verifications.
To verify and quantify our potential model, we performed density functional theory (DFT) calculation using CASTEP. The generalized 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 van der Waals force was also examed by adopting the DFT-D, its effect to the stability of the ice structures has reported, a small effect of red-shifts (~2 meV) for the main peaks of the vibrational spectra was observed. The phonon spectra were calculated using the CASTEP module with finite displacement method. The force constants produced by the CASTEP for the phonon spectra calculations were obtained from the output files. The force for the atom i, f i (= dE/dr i ), is simply the derivative of the total energy E. By applying a second small displacement for the atom j, the force constant k ij (= dU/dr i dr j ) is obtained for the pair atoms i and j. A force constant matrix k for the unit cell is a 3 N× 3 N matrix (where N is total number of atoms in the unit cell) and were constructed based on above procedure.
Ice VIII has a cubic structure containing two sets of interpenetrated ice sub-lattices as shown in Fig. 2a. x H shrinks and x C extends in stage I. As pressure increases, x H and x C both shrink in stage II like normal material.
The DFT calculation can reproduce phonon spectrum for Ice VIII as shown in Fig. 3a. The simulation reproduces the main features of inelastic neutron scattering (INS) spectrum 24 , such as the three peaks in the translational region (< 50 meV) and the two sharp peaks 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. Force constants k C , k H and k OO were also obtained from the calculation of phonon spectrum of the ice VIII and are plotted as a function of pressure in Fig. 3b. Apart from the trend for the k C , k B to merge at x OO    Potential Model Parametrization. Fitting from DFT-derived k-x results of stage I of P-dependent IceVIII (Fig. 3b), we get the relations among k i and x i (I = H, C, OO) according to our potential model, where we add negative sign on k H since DFT only take the positive square root value. Based on k-x curve, we fit the expression of f-x and V-x curves (Supplementary Eq. S1 and S2). The shapes of f-x and V-x are plotted in Fig. 4. In stage I, as f p increases, f H = f C should also increases according to Eq.(3). Thus, x H decreases while x C increases as indicated by black arrows in Fig. 4a. In stage II, when x H decreases to about 1.1 and reach to the curvature turning point, f H changes to decrease and x C will decreases back (as shown by blue arrows) to make f C = f H . The absolute value of repulsive force f OO always increases as distance x OO decreses. We focus on the parametrization of stage I in this work since it is more interesting than the normal compression-contraction in stage II.
Hence, for a given f P , there is only one set of solutions of x H and x C obtained in Fig. 5. Figure 5a shows that, for a given x H , f H = f C gets the x C . Then, Fig. 5b shows that, f OO can be calculated by x H + x C and f P is correlated with x H by f p = f OO − f H . Hence, we have obtained the f P -x H curve in Fig. 5b, from Figure 5b also indicates that, as f P increases, f OO rises much more significantly than f H . Thus, f OO is indeed a hidden force that should not be ignored.
Based on our potential model, foo, f H and fp are functions of x H and x C determined by Supplementary Eq. S1 and Fig. 4. We run the time-dependent dynamic process of external pressure on O:H-O bond in Fig. 6. The dynamic equation can be expressed as: = . The initial condition of x H and x C are the lengths (Å) without pressure. The partial differential equation is solved by 4 th order Runge-Kutta method (Supplementary Method). Figure 6 shows the x H and x C relaxation with the same initial lenghs and different external forces under f p = 0.1, 0.5 and 1.5 eV/Å. The initial lengths are bond lengths at fp = 0, which give energies to atoms to oscillate. If the process is slow, i.e. initial coordinates are around the stable point, the oscillation will decrease. x H and x C under fp = 0.5 and 1.5 eV/Å were taken average as ∑ x t t ( ())/