Structural Properties of Liquid Water and Ice Ih from Ab-initio Molecular Dynamics with a Non-local Correlation Functional

Equilibrium Born-Oppenheimer molecular dynamics simulations have been performed in the canonical ensemble to investigate the structural properties of liquid water and ice Ih (hexagonal ice) at 298 and 273 K, respectively, using a state-of-the-art non-local correlation functional, whilst size effects have been examined explicitly in the case of liquid water. This has led to improved agreement with experiments for pair distribution functions, in addition to molecular dipole moments, vis-à-vis previous flavours of ab-initio molecular dynamics simulation of water, highlighting the importance of appropriate dispersion. Intramolecular geometry has also been examined, in addition to hydrogen-bonding interactions; it was found that an improved description of dispersion via non-local correlation helps to reduce over-structuring associated with the Perdew-Becke-Ernzerhof (PBE) and other commonly-used functionals.


Introduction
The simulation of water [1][2][3][4][5][6][7] and ice [8], or indeed of clathrate hydrates [9], and other condensed matter aqueous systems in the liquid or solid state, using density functional theory-based functionals, can be quite problematic, owing in large part to challenges in modelling dispersion effectively [1][2][3][4][5][6][7][8][9].From the perspective of ab-initio molecular dynamics (AIMD), to achieve qualitatively similar results with experimental structural quantities, such as pairwise radial distribution functions (RDFs) for liquid water, it has been necessary to use temperatures perhaps 20% higher, due to over-structuring by the functionals [1], e.g., Perdew-Becke-Ernzerhof (PBE) [10].State-of-the-art simulation of ice Ih involving DFT and classical molecular dynamics has also been applied to compute RDFs quite accurately by Daskalakis and Hadjicharalambous [11], while Hernández de la Peña and Kusalik have studied quantum effects in both liquid water and ice using a variety of rigid water models, to extract RDFs [12].
However, recently, the first-principles van der Waals (vdW) density functional approach (vdW-DF) of Dion et al. has yielded impressive initial results in a more accurate and satisfactory treatment of dispersion, ipso facto [13], and this has proven promising in initial AIMD simulation of liquid water [14], particularly vis-à-vis addressing the problem of over-structuring associated with PBE and other commonly-used functionals with a less sophisticated treatment of dispersion.
In this study, we adopt the implementation of vdW-DF by Román-Pérez and Soler [15], using revPBE exchange in conjunction with non-local vdW correlation [13], termed "DRSLL" following relatively positive initial results in the simulation of liquid water [14], to carry out equilibrium Born-Oppenheimer molecular dynamics (BOMD) in the canonical ensemble of liquid water and ice Ih.Our goal is the study of structural properties of water and ice, particularly with respect to assessing the relative accuracy of DRSLL with respect to experiment.In assessing structural properties, we consider not only RDFs, but also the important matter of molecular dipole moment distributions (given the interplay between structural and induced polarisation features to shape dipole moments).Naturally, understanding system size effects on structural and dipolar properties is of relevance to ab initio simulation, in setting the "lower bounds" of what is necessary to capture physically accurate data from AIMD for liquid water and ice simulation.

Simulation Methodology
In this study, we used revPBE exchange in conjunction with non-local vdW correlation [13], within the pseudopotential localised basis set approximation, as implemented in the SIESTA code [16].Split valence plus polarisation (DZP) basis sets were employed for oxygen and hydrogen atoms, as in ref. [11].The Coulomb potential was expanded in a plane-wave basis with an energy cut-off of 300 Ryd; a reasonably large plane-wave expansion is important in order to ensure convergence in the calculated total energy and atomic forces.A Nosé-Hoover thermostat was applied with a period of 0.25 ps [17,18].Following 2 ps of relaxation during which the total energy and potential energy were found to converge, as well as intramolecular geometries, BOMD was carried out under periodic boundary conditions (PBC) for 5 ps at 298 and 273 K with a time step of 1 fs for liquid water and ice Ih, respectively.This obviates the somewhat less satisfactory practice of using slightly higher temperatures (~20%) to reduce problems of over-structuring for PBE [1] and other commonly-used functionals.Earlier validation in the micro-canonical ensemble indicated that a time step of 1 fs was satisfactory.
To gauge explicitly the size dependence of the computed parameters, three system sizes were used for liquid water at a density of 1 g/cm 3 , corresponding roughly to experimental ambient pressure conditions: 32, 60 and 116 molecules, all in a cubic simulation box (with respective dimensions of 9.864, 12.183, 15.142 Å), relaxed previously from classical MD [17] with the SPC/E model [19].Two orthorhombic unit cells containing 96 and 192 molecules of Hayward and Reimers [20] were used as starting configurations for hexagonal ice (Ih), in a proton-disordered structure, with a zero net dipole and a zero net quadrupole moment.The x, y and z box lengths were 13.521, 15.613 and 14.72 Å, respectively, for the 96-molecule system, with respective values of 18.028, 23.419 and 14.72 Å for the 192-molecule case.This gave a density very close to the experimental value of 0.92 g/cm 3 in both cases, at 273 K and ambient pressure [21].

Results and Discussion
For water, the pressure was found to be around −0.91 ± 0.57, −0.59 ± 0.49 and −0.48 ± 0.42 kbar for the 32-, 60-and 116-molecule systems, respectively.This is in general accord with the findings of Wang et al. [14], although the system size for AIMD used in that study comprised 64 molecules.In any event, this pressure is closer to the "desired" experimental value of ambient pressure than PBE [1][2][3][4][5][6], and this represents an advance, in some sense, over previous attempts at dispersion correction in water [22], which used empirical pair potentials proposed by Grimme [23] to account for vdW interactions, to essentially "recover" the pressure density relationship.For ice, the pressure was found to be −0.584± 0.486 kbar and −0.472 ± 0.424 kbar for the 96-and 192-molecule systems, respectively, which is relatively good in terms of (increasing) closeness to ambient pressure.This increasing closeness of the pressure vis-à-vis ambient pressure upon increasing system size for both ice and liquid water underlines the importance of achieving a "certain" threshold minimum size for the convergence of the pressure tensor (whatever the amplitude of the temporal fluctuations thereof, related to the inverse square root of the number of particles [17]).This important point was also highlighted by Wang et al. in their use of classical MD with the TIP4P (four-point transferable intermolecular potential) potential to illustrate the convergence of pressure estimates with liquid water systems of up to around a thousand molecules [14]; of the order of one hundred molecules was found necessary to achieve this.This is not unexpected, given that the corresponding cubic simulation box size is around 15 Å, and half of this would correspond to 7.5 Å, which is slightly larger than the "rule of thumb" of 2.2-times the Lennard-Jones 12-6 radius for a minimum cut-off (if any) for van der Waals interactions typical for MD under PBC; use of a box size appreciably smaller than around 15 Å would be expected to lead to truncations of "missing" pairwise dispersion interactions, especially in the attractive "basin" (which may be modelled empirically by, e.g., a Lennard-Jones form).This helps to explain the need for a system containing of the order of a hundred water molecules as being a reasonable, now computationally-feasible minimum size for AIMD to achieve basic convergence of pressure and energy terms in condensed matter aqueous systems (e.g., ice or liquid aqueous phases).
Given that the (admittedly tentative) rate of change in system pressure upon increasing system size declines for liquid water (albeit difficult to conclude for ice, given the use of two system sizes) and bearing in mind the above comments in relation to around a hundred water molecules being a reasonable "lower bound" for the system size, this would lead one to expect the physical results for the 116-molecule liquid water system to be the most reliable and that the 96-molecule ice Ih system may be somewhat reasonable in terms of the prediction of the structure.Below, however, we shall see that there are nuances and subtleties to this conclusion in the case of 96-versus 192-molecule ice.Bearing this in mind, the O-O RDFs are presented in Figures 1 and 2 for liquid water and ice, respectively, compared to experimental data [24,25], with results for the three different system sizes shown in Figure 1 in the case of water.Broadly, the water results in Figure 1 agree semi-quantitatively with careful calculations of the RDF from classical MD (with good quality empirical pairwise potentials) in liquid water in [12]; they are also in good accord with the DRSLL findings of Wang et al. [14], although Wang et al. also computed the case for DRSLL-PBE.Given the very good agreement of DRSLL in the present work with that of Wang et al. and Zhang et al. [26], we would also conclude that DRSLL-PBE is a reasonable model for structural properties of liquid water vis-à-vis DRSLL, albeit with over-correction for density and overly strong hydrogen bonds in the case of DRSLL-PBE [14].Del Ben carried out 64-molecule liquid water Monte Carlo simulations using MP2 [27], and the RDF in Figure 1 from DRSLL is in reasonable accord with MP2 results also (given MP2's propensity to model hydrogen bonds more accuracy vis-à-vis DFT, in general).It can be seen that the RDF in Figure 1 is in closer general agreement with the experimental data for the larger systems, particularly the 116-molecule case; Wang et al. were unable to simulate a system of this size, carrying out a 64-molecule ab initio simulation [14], while Zhang et al. [26] performed 32-molecule simulations.This somewhat closer agreement is especially evident beyond the immediate (hydrogen bonded) nearest neighbours; this comes as little surprise, given that the box length for the smallest 32-molecule system is only ~9.86 Å, meaning that obtainable pairwise interactions extend little beyond around 5 Å, hardly past the second coordination layer.It can be seen that DRSLL results in less structured liquid water compared to the experiment (and also vis-à-vis PBE [14], although that functional itself is over-structured).For ice, Figure 2 (for the 96-molecule system) shows that DRSLL leads also to a loss of structural features in comparison to the experiment, especially in terms of the secondary coordination shell and beyond.However, the ice RDFs are in good agreement with reference model calculations from both classical MD [11,12] and B3LYP [11].In Figure 2, the O-O RDF of the 192-molecule ice system is in near-identical agreement with the 96-molecule one, emphasising again the apparent point that around a hundred molecules suffice for a reasonable description of structure (at least from an RDF standpoint) in condensed-matter aqueous phases.The relatively small negative pressures of both liquid water and ice for the given experimental densities (the magnitude of which declines, reassuringly, with increasing size) suggest that refinements to the revPBE local exchange may help in correcting this, so as to get better agreement for the position of the first peak.The peak heights are in reasonable accord with experiment in Figures 1 and 2 for liquid water and ice, although the width of the histogramming used in the computation of the RDF can of course affect this to some extent, making exact quantitative comparison somewhat more difficult.DRSLL can be seen to over-structure the RDF of liquid water in Figure 1, whilst it essentially under-structures it for ice Ih in Figure 2; given the non-empirical nature of the DRSLL model, it has not been parameterised specifically for the liquid or ice state.As such, under-or over-structuring of either phase is possible (and indeed observed in this case).
To give a reference point for the effect of the liquid/ice Ih environment, the O-H bond length and the H-O-H bond angle for an isolated water monomer were 0.9714 Å and 104.53°, in good accord with [26].Furthermore, the O-O length for a hydrogen-bonded dimer was found to be 3.082 Å, in good agreement with [14].
In terms of intramolecular structure (using the 116-molecule system for water), it was found that the O-H bond lengths were 0.978 ± 0.012 and 0.983 ± 0.008 Å for water and (96-molecule) ice, whilst the respective H-O-H angles were 106.3° ± 1.4 and 108.7° ± 0.8.There were 3.77 ± 0.08 hydrogen bonds per water molecule in the case of the liquid, with hydrogen bond lengths of 2.00 ± 0.14 Å and 1.76 ± 0.04 Å in water and ice, respectively.This is in general accord with experimental results for hydrogen bond lengths and intramolecular geometry in liquid water and ice [21].It would appear that DRSLL models this aspect of water structure quite well, with the main problem being under-structuring in the second coordination shell.There was relatively good agreement for intramolecular geometry for liquid water with respect to system size, whereas the bond angle and O-H bond length were both ~0.5% larger for 192-molecule ice vis-à-vis the smaller system, indicating that while intermolecular geometry (as characterised by the O-O RDF in Figure 2) is relatively insensitive to doubling of ice system size, there are more subtle intramolecular effects, not as evident in the case of liquid water.Given the crystalline structure of ice Ih, as opposed to the lack of longer range order in liquids, this tends to reinforce longer range local electronic polarisation and geometry effects.Given this observed more subtle effect in ice upon "scale-up" that is not as evident in the liquid state, the dipole moment distribution was computed from AIMD, in an effort to evince and clarify putative system size effects.This was done by the Hirshfeld scheme for partial charge decomposition [28] (of AIMD trajectories).The monomer dipole moment was found to be 1.80 D, in good agreement with the experimental value of 1.85 D [29].For reference, experimental values of the molecular dipole moment in liquid water are 2.95 ± 0.6 D (inferred from X-ray data) [30], whilst in ice Ih, it is ~3.09D (inferred from an induction model approach) [31].In Figure 3, the most probable value of the dipole moment has "converged" to ~2.85 D for 116-molecule liquid water, with an upward shift of some ~0.12 D from 2.73 D for the 32-molecule case (with a similar spread).The intermediate, 64-molecule distribution is centred at ~2.83 D, very close to the 116-molecule case, suggesting that the 32-molecule system is lacking the longer range electrostatic effects needed to induce greater polarisation in water characteristic of the bulk state.This DRSLL value of ~2.85 D is in reasonable accord with a DRSLL-PBE value of 2.93 ± 0.012 D for liquid water [26] and the experimental value of 2.95 ± 0.6 D [30].In contrast, Figure 4 depicts a more subtle effect in ice, with both system sizes having essentially the same most probable value of ~3.18 D (close to the experimental mean value of ~3.09 D [31]), but the larger, 192-molecule system showing evidence of a more narrow spread and distribution.The greater system size and extent of "real" (as opposed to "replica") crystalline periodicity appears to "reinforce" induced polarisation in a more subtle way than the smaller system, restricting the spread of the distribution.Indeed, the typical behaviour of decline in magnitude of the standard deviation (or "fluctuation") with the square root of the system size [17] may serve to explain this in the case of ice, for which the smaller, 96-molecule-sized system has already essentially converged the centre point of its distribution towards ~3.18 D. However, in the liquid case (cf. Figure 3), the spreads are similar, as the system is still reaching the centre point of the distribution as it increases in size.

Conclusions
Equilibrium BOMD simulation of liquid water and ice in the canonical ensemble has been performed to investigate some of the salient structural properties of liquid water and ice, with a view toward assessing an improved description of the dispersion via a non-local correlation with the DRSLL functional and establishing a "lower bound" on reasonable system size.It was found that DRSLL offers an improvement over PBE and other functionals, in the sense that the intended physical temperature can be used and over-structuring is avoided, as is often the case with PBE and other commonly-used functionals, which do not include dispersion interactions in what would generally be regarded a satisfactory manner; these mirror largely the findings of Wang et al. [14] in the case of liquid water.However, a reasonable "threshold" size of around a hundred molecules is desired for AIMD under PBC, to allow the pressure, energy and molecular dipole moment distribution to converge and offer more adequate agreement with the experiment; the increasing availability of linear scaling DFT (e.g., in SIESTA or other implementations [32]) means that AIMD is increasingly able to probe systems containing several hundred atoms systematically.However, for this to be meaningful, rigorous incorporation of dispersion interactions is a key desideratum, particularly for aqueous systems.Although DRSLL has offered a certain level of improvement, refinements to the revPBE local exchange may well offer the scope for future improvements, as would developments in the underlying approaches to vdW-DF [33].

Figure 1 .Figure 2 .
Figure 1.O-O radial distribution function (RDF) for liquid water at 298 K.The ambient pressure experimental data is from[24], whilst results for the three different system sizes are presented for DRSLL.

Figure 3 .
Figure 3. Probability distribution of the molecular dipole moment for DRSLL liquid water at 298 K, for the smallest (32-molecule) system on the left and the largest (116-molecule) on the right.

Figure 4 .
Figure 4. Probability distribution of the molecular dipole moment for DRSLL ice Ih at 273 K, with the smaller (96-molecule) system "below" exhibiting a wider distribution and the more narrow distribution for the large one (192-molecule) evident "above".