On the Accuracy of van der Waals Inclusive Density-Functional Theory Exchange-Correlation Functionals for Ice at Ambient and High Pressures

On the Accuracy of van der Waals Inclusive Density-Functional Theory Exchange-Correlation Functionals for Ice at Ambient and High Pressures

Biswajit Santra    Jiří Klimeš    Alexandre Tkatchenko    Dario Alfè    Ben Slater    Angelos Michaelides angelos.michaelides@ucl.ac.uk    Roberto Car    Matthias Scheffler Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, 14195 Berlin, Germany
Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA
Faculty of Physics and Center for Computational Materials Science, University of Vienna, Sensengasse 8/12, A-1090 Wien, Austria
London Centre for Nanotechnology, University College London, London WC1E6BT, UK
Department of Physics and Astronomy, University College London, London WC1E6BT, UK
Department of Earth Sciences, University College London, London WC1E6BT, UK
Department of Chemistry, University College London, London WC1E6BT, UK

Density-functional theory (DFT) has been widely used to study water and ice for at least 20 years. However, the reliability of different DFT exchange-correlation (xc) functionals for water remains a matter of considerable debate. This is particularly true in light of the recent development of DFT based methods that account for van der Waals (vdW) dispersion forces. Here, we report a detailed study with several xc functionals (semi-local, hybrid, and vdW inclusive approaches) on ice I and six proton ordered phases of ice. Consistent with our previous study [Phys. Rev. Lett. 107, 185701 (2011)] which showed that vdW forces become increasingly important at high pressures, we find here that all vdW inclusive methods considered improve the relative energies and transition pressures of the high-pressure ice phases compared to those obtained with semi-local or hybrid xc functionals. However, we also find that significant discrepancies between experiment and the vdW inclusive approaches remain in the cohesive properties of the various phases, causing certain phases to be absent from the phase diagram. Therefore, room for improvement in the description of water at ambient and high pressures remains and we suggest that because of the stern test the high pressure ice phases pose they should be used in future benchmark studies of simulation methods for water.

I Introduction

Density-functional theory (DFT) is now widely used to study water and ice in a range of different environments, including for example bulk water, water at interfaces, and water under confinement. Most DFT studies of water have involved the application of semi-local generalized gradient approximations (GGA) for the exchange and correlation (xc) energy. Whilst these studies have proved to be very useful in providing insights into the structure and properties of water, there are persistent question marks over the quantitative accuracy of such xc functionals, in particular for the treatment of condensed phase water which is held together by hydrogen (H) bonding and van der Waals (vdW) interactions. Over the years this has prompted a number of benchmark studies focused on gas phase water clusters Santra et al. (2007, 2008, 2009); Dahlke et al. (2008); Xu and Goddard III (2004); Su et al. (2004); Xantheas (1995); Kim and Jordan (1994); Anderson and Tschumper (2006); Shields and Kirschner (2008); Csonka et al. (2005); Dahlke and Truhlar (2005); Gillan et al. (2012); Ireta et al. (2004); Tsuzuki and Lüthi (2001); Novoa and Sosa (1995); Hammond et al. (2009); Wang et al. (2010), liquid water, Grossman et al. (2004); Schwegler et al. (2004); Fernández-Serra and Artacho (2004); McGrath et al. (2005); McGrath et al. (2006a, b); Lee and Tuckerman (2006, 2007); Todorova et al. (2006); VandeVondele et al. (2005); Guidon et al. (2008); Kühne et al. (2009); Sit and Marzari (2005); Fernández-Serra et al. (2005); Asthagiri et al. (2003); Mattson and Mattson (2009); Kuo et al. (2004); Sprik et al. (1996); Silvestrelli and Parrinello (1999); Yoo et al. (2009); Zhang et al. (2011a); Chen et al. (2003); Morrone and Car (2008); Laasonen et al. (1993); Akin-Ojo and Wang (2011); Alfè et al. (2013) and crystalline phases of ice. Santra et al. (2011); Murray and Galli (2012); Hamada (2010); Kolb and Thonhauser (2011); Kambara et al. (2012); Fang et al. (2013); Gillan et al. (2013); Pamuk et al. (2012); Lee et al. (1992); Hamann (1997); Singer et al. (2005); Tribello et al. (2006); de Koning et al. (2006); Feibelman (2008); Hermann and Schwerdtfeger (2008); Pan et al. (2008); Erba et al. (2009); Morrone et al. (2009); Militzer and Wilson (2010); Labat et al. (2011)) While we currently have a relatively clear understanding about the performance of various  xc functionals for gas phase clusters, this is far from being established for ice and liquid water. This is particularly true in light of recent work which has shown that vdW dispersion forces are important for the accurate description of different properties of water. Santra et al. (2011); Murray and Galli (2012); Hamada (2010); Kolb and Thonhauser (2011); Kambara et al. (2012); Fang et al. (2013); Gillan et al. (2013); Lin et al. (2009); Schmidt et al. (2009); Wang et al. (2011); Møgelhøj et al. (2011); Zhang et al. (2011b); Jonchiere et al. (2011); Yoo and Xantheas (2012); Ma et al. (2012); zha ()

Understanding the role of vdW forces in water has been greatly helped by the emergence of various approaches for accounting for vdW forces within the framework of DFT. Dion et al. (2004); Tkatchenko and Scheffler (2009); Klimeš et al. (2010); von Lilienfeld et al. (2004); Silvestrelli (2008); Sato and Nakai (2009); Grimme et al. (2010); Vydrov and Voorhis (2009); Becke and Johnson (2005); Tao et al. (2012) In the last few years many of the vdW inclusive DFT xc functionals have been used to investigate the effects of vdW on the structural, energetic, and vibrational properties of liquid water. Lin et al. (2009); Schmidt et al. (2009); Wang et al. (2011); Møgelhøj et al. (2011); Zhang et al. (2011b); Jonchiere et al. (2011); Yoo and Xantheas (2012); Ma et al. (2012); zha () Overall, with vdW inclusive xc functionals there are indeed improvements in certain calculated properties of liquid water. For example, the first peak in the oxygen-oxygen radial distribution function is generally reduced and brought into closer agreement with experiment. However, the accuracy of the computed properties strongly depends on the methods chosen to incorporate vdW as well as the technical details of the molecular dynamics simulations. There is, of course, also the challenge of accounting for quantum nuclear effects, which is rarely done in ab initio studies of liquid water. Chen et al. (2003); Morrone and Car (2008) However, in contrast to liquid water, the various crystalline phases of ice represent a relatively straightforward set of structures against which DFT methods can be tested. Indeed there are at present 15 experimentally characterized ice phases with water molecules in a number of distinct arrangements, H bond networks, and densities. Petrenko and Whitworth (2003); Salzmann et al. (2009); sal (); Salzmann et al. (2011) Many of the ice phases are complicated with disordered water arrangements (so called “proton disordered”). However, some phases have relatively simple proton ordered arrangements of water molecules, and it is these phases that are particularly suitable as benchmarks. Furthermore, thanks to Whalley’s extrapolations of the experimental finite temperature and pressure phase coexistence lines to zero temperature, for some of these phases there are even estimates of the internal energy differences, Whalley (1984) which makes theoretical benchmarks particularly straightforward and mitigates the need for expensive free energy calculations.

Figure 1: Unit cells of the ice phases (I, IX, II, XIII, XIV, XV, and VIII) studied here. The ice I structure (96 water molecule) is obtained from ref. Pan et al. (2008) and all the proton ordered phases of ice are obtained from various scattering experiments. Line and Whitworth (1996); Londono et al. (1993); Lobban et al. (2002); Salzmann et al. (2006, 2009); Kuhs et al. (1984) The optimized coordinates of the ice structures are given in the supplementary material. sup ()

In an earlier study on ice we found that the effects of vdW become increasingly important upon going from the low- to the high-density phases and capturing this variation in the vdW energy is essential to get the transition pressures between the ice phases within an order of magnitude of experiment. Santra et al. (2011) Here, we extend on the previous study significantly by reporting results on the accuracy of the cohesive properties of individual phases of ice obtained from a wide range of vdW inclusive functionals. Also by looking at the enthalpies of ice as a function of pressure we have obtained a more detailed picture of the stability range of each ice phase predicted from the different functionals. The approaches used here include: i) vdW, which involves an explicit summation of pair-wise (two-body) vdW dispersion interactions among all atom pairs using their respective vdW coefficients which are functionals of the electron density; Tkatchenko and Scheffler (2009) (ii) vdW, an extension of vdW that accounts for electrodynamic screening and many-body vdW interactions within the dipole approximation; Tkatchenko et al. (2012) and (iii) various functionals from the “vdW-DF” family. Dion et al. (2004); Lee et al. (2010); Klimeš et al. (2010) All vdW-DFs are calculated via a model dynamic response function and long range pairwise approximations. Dion et al. (2004) These various vdW inclusive approaches have been reasonably successful in modeling a wide variety of materials Langreth et al. (2009); Tkatchenko et al. (2010); Klimeš et al. (2011); DiStasio Jr. et al. (2012); Klimeš and Michaelides (2012); Reilly and Tkatchenko (2013) including different phases of water starting from clusters Santra et al. (2008); Kelkkanen et al. (2009); Klimeš et al. (2010) to condensed phases. Santra et al. (2011); Murray and Galli (2012); Wang et al. (2011); Møgelhøj et al. (2011); Zhang et al. (2011b); zha () In this study we find that all vdW inclusive functionals considered predict phase transition pressures in much better agreement with experiment than the functionals which do not include vdW. However, the precise values of the lattice constants and lattice energies are highly sensitive to the choice of vdW inclusive method. Moreover, none of the functionals can simultaneously produce energetics and volumes of the ice phases with high enough precision to yield a phase diagram that correctly captures all the phases found in experiments.

In the next section details of the simulation methods are provided. This is followed by discussions of our results for the equilibrium lattice energies (III.1), the equilibrium volumes (III.2), enthalpies (III.3), and a decomposition of total energies focusing on exchange and correlation energies (III.4). Conclusions and a short perspective on future work are given in section IV.

Ii Simulation Details

We have computed and analyzed the equilibrium lattice energies, volumes, and enthalpies of several ice phases. This includes the ambient pressure phase of ice, ice I, and all the proton ordered high-pressure phases, namely, in order of increasing pressure, ice IX, II, XIII, XIV, XV, VIII. We have focused on proton ordered phases because they are more straightforward to model than the proton disordered phases. The initial structures used for the proton ordered phases have been obtained from experiment Line and Whitworth (1996); Londono et al. (1993); Lobban et al. (2002); Salzmann et al. (2006, 2009); Kuhs et al. (1984) and the unit cells used are shown in Fig. 1. Proton disordered ice I is modeled with the 12 water unit cell proposed by Hamann. Hamann (1997) The results obtained from the 12 molecule cell have also been compared to results from a unit cell of 96 water molecules (Ref. Pan et al. (2008, 2010)). These results reveal that the lattice energies obtained from the 12 and 96 water molecule unit cells are within 1 meV/HO and the equilibrium volumes differ only by 0.01 Å/HO with PBE. Perdew et al. (1996)

The lattice energy per HO () of ice is obtained by subtracting the total energy of isolated HO molecules () from the total energy of the ice unit cell () containing molecules of HO, i.e.,


At zero pressure the theoretical equilibrium lattice energies and volumes are obtained by varying the lattice parameters isotropically within 20% of the experimental values and fitting the resultant energy-volume curves to the Murnaghan equation of state. mur () By isotropic variation we mean that the ratios of the lattice parameters are kept fixed at the experimental value, which is a reasonable approximation that has an insignificant influence on the computed properties. For example, performing a rigorous test on ice VIII by varying the ratio of the lattice parameters provides changes of 0.5 meV/HO and 0.02 Å/HO, respectively in the equilibrium lattice energy and volume when compared to the results obtained by fixing the ratio at the experimental value. san () Also previously it was shown that for ice I the equilibrium ratio is very similar (within 0.4%) to the experimental value when calculated with various xc functionals. Feibelman (2008)

The properties of the various ice phases have been investigated with seven functionals, representing a number of different classes of functional. These include, PBE, a widely used GGA functional, and PBE0, Adamo and Barone (1999) a hybrid exchange variant of PBE. Neither of these functionals account for vdW forces. We have also considered PBE+vdW and PBE0+vdW, vdW inclusive versions of PBE and PBE0, where the vdW interaction is calculated with the Tkatchenko and Scheffler (TS) scheme. Tkatchenko and Scheffler (2009) The xc energy () in this scheme takes the form


where, is the PBE or PBE0 exchange, is the LDA correlation, is the PBE semi-local correlation correction, and is the vdW energy in the TS scheme. In addition, we have employed an extension of the approach which takes into account many-body (MB) dispersion and long-range electrostatic screening. Tkatchenko et al. (2012) In this case in Eq. (1) is replaced with MB dispersion energy terms and the xc functional is referred to as PBE0+vdW. Another approach to incorporate vdW within DFT is employed here, specifically the approach generally referred to as “vdW-DF”. Dion et al. (2004) In this case the total xc energy takes the form


where, is GGA exchange, and is the nonlocal correlation energy through which the vdW interactions are captured. We have used three functionals from this category which we refer to as revPBE-vdW, Dion et al. (2004) optPBE-vdW, Klimeš et al. (2010) and rPW86-vdW2. Lee et al. (2010) The difference between revPBE-vdW (originally proposed in ref. Dion et al. (2004)) and optPBE-vdW is in the exchange functional only. The former employs revPBE Zhang and Yang (1998) exchange, whereas the latter uses optPBE exchange Klimeš et al. (2010) which was developed by fitting interaction energies obtained for the S22 data set. Jurečka et al. (2006) Typically optPBE exchange is less repulsive than revPBE at intermediate and short inter-atomic distances. Compared to the above two functionals rPW86-vdW2 has a different exchange rPW86 Murray et al. (2009) and a modified nonlocal correlation functional. Lee et al. (2010) We should note that all three vdW-DF functionals utilize GGA exchange.

The calculations with PBE and PBE+vdW were performed with the all electron numeric atom-centered orbital (NAO) basis set code FHI-aims. Blum et al. (2009) Sufficiently large basis sets (“tier2” for H and “tier3” for O) were employed to calculate total energies and to optimize structures. PBE0, revPBE-vdW, optPBE-vdW, and rPW86-vdW2 calculations were done with the VASP code Kresse and Hafner (1993); Kresse and Furthmüller (1996) with the hardest projector-augmented wave (PAW) pseudopotentials and a 1000 eV plane-wave basis set cut off. The is calculated self-consistently with the efficient algorithm of Román-Pérez and Soler Román-Pérez and Soler (2009) employing 30 interpolation points for the function with a saturation value a.u., as implemented by Klimeš et al. in VASP. Klimeš et al. (2011) These settings are found to be very accurate and details of the implementation and tests performed for a variety of solids can be found in ref. Klimeš et al., 2011. For all the ice structures the atoms are fully relaxed with all of the xc functionals (except with PBE0+vdW and PBE0+vdW) without any symmetry constraints until all forces are less than 0.01 eV/Å. The energy-volume curves of each ice phase with PBE0+vdW and PBE0+vdW were produced by performing single point energy calculations on the PBE0 optimized geometries at different volumes. For the calculations of any GGA exchange based functional the number of k points are chosen so that the spacing in the k point grid in each direction of reciprocal space is within 0.02 Å to 0.04 Å for all of the ice phases. For the hybrid functional (PBE0) calculations the number of k points are doubled in each direction compared to the GGA calculations, which provides total energies converged to within 1 meV/HO. With VASP the energy of the water monomer was calculated within a cubic cell of length 20 Å.

Iii Results

Figure 2: Differences in the lattice energies of ice calculated () with various DFT functionals and DMC Santra et al. (2011) compared to experiment (). Zero on the Y axis indicates perfect agreement with experiment.
Expt. -610 (0) -606 (5) -609 (1) -577 (33)
DMC -605 (0) -609 (-4) -575 (30)
PBE -636 (0) -587 (49) -567 (69) -556 (80) -543 (93) -526 (110) -459 (177)
PBE0 -598 (0) -557 (41) -543 (55) -530 (67) -518 (80) -504 (94) -450 (148)
PBE+vdW -714 (0) -705 (9) -698 (16) -695 (19) -690 (24) -678 (36) -619 (95)
PBE0+vdW -672 (0) -670 (2) -666 (6) -661 (11) -656 (16) -646 (26) -596 (76)
PBE0+vdW -672 (0) -663 (9) -656 (16) -648 (22) -642 (30) -629 (43) -589 (83)
revPBE-vdW -559 (0) -563 (-4) -556 (3) -555 (4) -552 (7) -545 (14) -517 (42)
optPBE-vdW -668 (0) -673 (-5) -667 (1) -666 (2) -664 (4) -656 (12) -630 (38)
rPW86-vdW2 -619 (0) -621 (-2) -618 (1) -615 (4) -605 (14) -605 (14) -586 (33)
Ref. Whalley, 1984; The DMC statistical error bar is 5 meV/HO (Ref. Santra et al., 2011)
Table 1: Equilibrium lattice energies of different ice phases with various methods. The relative lattice energies of the high-pressure ice phases with respect to ice I are given in parenthesis. All energies are in meV/HO.

In this section we report how the above mentioned DFT xc functionals describe the different phases of ice by examining equilibrium lattice energies, volumes at zero pressure and the relative enthalpies of the various phases. Subsequently we report an analysis of the individual contributions from exchange and correlation to the lattice energies.

iii.1 Lattice energies at zero pressure

The lattice energy is one of the key characteristic quantities of a solid and we use it here to evaluate the performance and deficiencies of different xc functionals in describing ice. Previously most analysis of lattice energies concentrated on ice IHamann (1997); Feibelman (2008); Hamada (2010) however, recently we showed that studying ice I alone is not enough to establish the general behavior of an xc functional over the entire phase diagram of water. Santra et al. (2011) Here we have calculated the lattice energies of different ice phases using a wide variety of xc functionals and made comparisons with experiments Whalley (1984); Röttger et al. (1994); Londono et al. (1993); Salzmann et al. (2006, 2009) and diffusion quantum Monte Carlo (DMC) Santra et al. (2011) whenever possible (Table 1). We note that the DFT and DMC lattice energies reported in Table 1 do not include nuclear zero-point energies (ZPEs) and are directly comparable with the experimental lattice energies provided by Whalley, Whalley (1984) in which ZPE contributions were removed and the energies were extrapolated to 0 K. In Fig. 2 the differences in the calculated and experimental lattice energies are shown for ice I, IX, II, and VIII for all the functionals considered. It can be seen that for the phases for which DMC data is available the agreement between DMC and experiment is excellent, differing only by 5 meV/HO at most, which is also the size of the DMC statistical errors.

Figure 3: Percentage differences in the calculated equilibrium volumes (V) compared to experiment (V) (a) without and (b) with zero point energies (ZPE) in V. The zero value on the Y axis designates the experimental reference values. Here positive errors indicate larger computed volumes and negative errors indicate smaller computed volumes compared to experiment.

We begin with the performance of functionals which do not account for vdW (PBE and PBE0) on the absolute values of the lattice energies. The behavior of PBE for ice I is well known, it over-estimates the lattice energy: overestimations of between 30 to 100 meV/HO have been reported depending on the computational set-up used (mainly the quality of the basis sets and pseudopotentials). Pan et al. (2008); Hamada (2010); Carrasco et al. (2011); Murray and Galli (2012) Here, using a full potential all-electron approach and very tightly converged NAO basis sets, we find that PBE over-estimates the lattice energy of ice I by 26 meV/HO (Fig. 2). This is in close agreement with the results from highly converged PAW calculations reported in ref. Feibelman, 2008. Interestingly the established notion that PBE overestimates the lattice energies of ice does not hold for the high-density phases. For example, PBE exhibits a 125 meV/HO underestimation for the lattice energy of ice VIII. The story is somewhat similar for PBE0, the hybrid variant of PBE. PBE0 predicts a very good lattice energy for ice I (only 15 meV/HO less than experiment) but simultaneously underestimates the lattice energy of the high-density ice phases. We believe that the behavior seen here for PBE and PBE0 is likely to apply to many other GGA and hybrid xc functionals. For example, our calculations show that BLYP and revPBE GGA functionals under-estimate the lattice energy of ice VIII by 246 meV/HO and 316 meV/HO, respectively. Similar findings have also been reported for B3LYP. Kambara et al. (2012)

In general we find that when vdW is accounted for the differences between the calculated and experimental lattice energies are much less sensitive to the particular phases being examined (Fig. 2). Considering first the vdW scheme, with the PBE+vdW and PBE0+vdW functionals the lattice energies are on average 100 meV/HO and 60 meV/HO, respectively too large compared to experiment. The smaller errors obtained from PBE0+vdW largely arise from the difference between PBE and PBE0 since the vdW contributions from these functionals obtained with the TS scheme are similar (to within 15 meV/HO). The contributions from vdW interactions beyond two-body vdW to lattice energies are found to be small for the phases considered. Specifically, the PBE0+vdW approach reduces the error by 7-17 meV/HO for the high-density phases compared to the standard PBE0+vdW. A noticeable exception to the consistent performance of vdW found for ice I, IX, and II is the highest density ice VIII phase (by 40-50 meV/HO [Fig. 2]). This inconsistency is largely due to the shortcomings of the damping function used in the vdW approach in describing the interpenetrating H bond network in ice VIII which has water molecules that do not form H bonds with each other as close as 2.9 Å apart.

When the vdW-DF functionals are utilized the errors are more consistent for the phases (Fig. 2). However, the magnitude and sign of the error varies considerably from one functional to another, e.g., on average optPBE-vdW produces too large (60 meV/HO) and revPBE-vdW produces too small (50 meV/HO) lattice energies compared to experiment. The fact that revPBE-vdW underestimates the lattice energy is not a surprise and consistent with results obtained with this functional for small molecules and water clusters. Gulans et al. (2009); Kelkkanen et al. (2009); Klimeš et al. (2010); Lee et al. (2010); Klimeš and Michaelides (2012) Previously reported lattice energies of ice I with revPBE-vdW are 30-35 meV/HO larger than what we obtain here. This is not a very substantial difference and we suspect it is mostly down to differences in pseudopotentials. Hamada (2010); Murray and Galli (2012) Here, rPW86-vdW2 provides the best agreement with experimental lattice energies being consistently within 15 meV/HO of experiment for all ice phases.

Without ZPE With ZPE
Expt. 32.05 25.63 24.97 23.91 23.12 22.53 20.09
DMC 31.69 24.70 19.46
PBE 30.79 26.11 25.01 24.08 23.27 22.82 20.74 1.69 0.57 4.00 2.99
PBE0 30.98 26.06 24.84 23.94 23.07 22.62 20.27 1.03 -0.14 3.14 2.32
PBE+vdW 29.67 23.86 23.62 22.44 21.71 21.47 20.13 5.52 -5.47 4.51 -3.49
PBE0+vdW 29.88 23.85 23.63 22.47 21.74 21.45 19.70 5.39 -5.39 4.05 -3.05
PBE0+vdW 29.42 23.87 23.26 22.26 21.45 21.10 18.90 6.88 -6.88 5.08 -5.08
revPBE-vdW 34.38 27.94 27.62 26.38 25.54 25.10 22.96 10.27 10.27 15.09 15.09
optPBE-vdW 31.63 25.50 25.15 23.99 23.20 22.75 20.55 0.92 0.42 3.21 3.09
rPW86-vdW2 33.69 26.65 26.35 25.07 24.24 23.74 21.27 5.09 5.09 8.22 8.22
10 K, Ref. Röttger et al., 1994; 30 K Ref. Londono et al., 1993; 0 K Ref. Whalley, 1984; 80 K Ref. Salzmann et al., 2006; 80 K Ref. Salzmann et al., 2009;
The DMC statistical errors are 0.01, 0.20, 0.02 Å/HO, respectively for ice I, II, VIII
(Ref. Santra et al., 2011); See Ref. not,a
Table 2: Comparisons of the calculated and experimental equilibrium volumes (Å/HO) of the various ice phases. MAE is mean absolute error (%) and ME is mean error (%) (averaged over all the ice phases) with respect to the experimental volumes. Errors with and without zero point vibration (ZPE) are shown. For the MAEs the positive sign indicates larger volumes and the negative sign smaller volumes compared to experiment.

The role played by vdW interactions in the phase diagram of ice is most evident when relative energies between the ice phases are considered. Both the vdW and vdW-DF approaches provide results which are in much closer agreement with experiment and DMC than the standard GGA and hybrid functionals in this regard. Table 1 shows that experimentally the energy difference between ice I and the highest density ice VIII phase is 33 meV/HO. Although the relative stabilities of ice XIII, XIV, and XV are not known (either from experiment or DMC) they should also fall within the 33 meV/HO window since ice VIII is the least stable phase at zero pressure and 0 K (of all the phases studied here). However, when calculated with PBE and PBE0 the energy difference between ice I and ice VIII is far too large (140 meV/HO). Likewise, when comparison with experiment is possible the phases between ice I and ice VIII are also destabilized too much. All vdW inclusive functionals reduce the energy differences between the phases, bringing them into much closer agreement with experiment. For example, the energy difference between ice I and ice VIII comes down to only 76 meV/HO and 33 meV/HO with PBE0+vdW and rPW86-vdW2, respectively.

As noted earlier, ZPE contributions were not considered in the above discussions as ZPE effects do not play a significant role in determining the relative energies of the various ice phases. The main effect from ZPE contributions is to reduce the lattice energies by about 120-110 meV/HO. This applies across the board for all phases and functionals considered, although a small monotonic decrease on the level of 6-10 meV/HO is seen upon going from the low- to the high-density phases. However, as we will see in the next section ZPE effects do influence the equilibrium volumes significantly.

iii.2 Equilibrium volumes

The equilibrium volume, which is a measure of the density of the phases, is another important quantity used to assess the performance of different theoretical methods. Previous efforts focused mostly on calculating the density of ice I, whereas here we seek to understand how functionals perform for a range of phases. Table II and Fig. 3 show comparisons of the calculated and experimental volumes. The equilibrium volumes of ice I obtained using both PBE and PBE0 are 4% smaller than experiment, which is consistent with previous calculations. Feibelman (2008); Murray and Galli (2012) Upon going to higher densities, however, contrasting behavior is observed and as one moves to higher densities there is an increasing tendency to overestimate the volume. Indeed for ice VIII the volume is overestimated by 4% with both PBE and PBE0. Clearly, the greater underestimation of vdW interactions at higher densities leads to a progressive overestimation of the equilibrium volumes. Overall though, and in contrast to the lattice energies, the performance of PBE and PBE0 for the equilibrium volumes are reasonable, differing by 2% from experiment when averaged over all ice phases.

Figure 4: Change in the enthalpies (H) of the ice phases relative to the enthalpy of ice I at T=0 and P=0 calculated with (a) PBE0, (b) PBE0+vdW, (c) PBE0+vdW, (d) revPBE-vdW, (e) optPBE-vdW, and (f) rPW86-vdW2. The vertical dotted lines indicate the transition pressures. The most stable ice phases along the pressure axis are indicated at the bottom of each panel. The insets show elaborations of the PBE0+vdW and optPBE-vdW plots within the 0.0-0.5 GPA pressure range.

With vdW the volumes are decreased by 3-9% from their parent functionals (PBE and PBE0) and in comparison to experiment the volumes actually become worse. The underestimated volumes obviously correlate with the underestimated lattice energies predicted by these approaches (Fig. 2). Going beyond two-body vdW the equilibrium volumes are reduced further and compared to experiment the average difference becomes 7%. It is noteworthy that vdW reduces the volume of ice VIII more than vdW, and as a result it improves the relative change in the volume with respect to ice XV, which in turn substantially affects the phase transition pressures (section III.3). Contrary to the performance of the vdW approaches, the equilibrium volumes obtained from revPBE-vdW are on average 10% too large (Fig. 3). Such behavior of revPBE-vdW has been attributed to the overly repulsive revPBE exchange functional and is analogous to what has been found earlier with this functional for many other solids. Klimeš et al. (2011); Klimeš and Michaelides (2012) The rPW86-vdW2 functional incorporates improvements in both the exchange and correlation components of the functional compared to revPBE-vdW and we find that this is reflected in improved volumes, being on average 5% larger than experiment. However, this performance is still inferior to PBE and PBE0. Of the vdW inclusive functionals optPBE-vdW provides the smallest average error being within 2% of experiment for all the ice phases.

Unlike the lattice energies, the effects of ZPE on the equilibrium volumes cannot be ignored, especially in the higher density phases. The effect of ZPE on the equilibrium volumes is estimated by computing free energy as a function of volume as , where being the frequency of phonon mode at a given volume. In line with the recent study of Murray and Galli, Murray and Galli (2012) we find that with PBE the volumes of ice I and VIII increase by 0.5% and 5.5%, respectively when ZPE effects are accounted for (Fig. 3(b) and Table 2). Indeed overall we find that the ZPE effects gradually increase from the low- to high-density phases and depending on the functional the increase in the equilibrium volume for the highest density ice VIII phase is somewhere between 3-6%. not (a) Thus, compared to experiments the mean absolute error in predicting volumes of phases increases by 2% for vdW-DF and decreases by 1.5% for vdW when ZPE effects are accounted for. Overall for all the vdW inclusive functionals, when ZPE effects are taken into account optPBE-vdW and PBE0+vdW are the two best functionals in terms of predicting volumes (Table 2).

iii.3 Enthalpy

Apart from absolute lattice energies and densities, accurate predictions of phase transitions are important if a functional is to be of real value in exploring the phase diagram of water. In ref. Santra et al. (2011) we showed that vdW interactions had a huge impact on the predicted phase transition pressures between the various phases of ice considered. Here we extend this study by calculating the enthalpies of different phases to establish the most stable phases at different pressures predicted by the various xc functionals. Pressures, P(V), at different volumes have been calculated from the Murnaghan equation of state: mur ()


where , , are the equilibrium bulk modulus, the derivative of the bulk modulus with respect to pressure, and the equilibrium volume at zero pressure, respectively. Fig. 4 shows the enthalpies of the various phases as a function of pressure relative to the enthalpy of ice I (at 0 K and zero pressure). The most stable phase at each pressure is the one with the lowest enthalpy and the crossovers between different phases indicate the pressures at which the phase transitions are predicted to occur. The phase transition pressures predicted from all functionals are also summarized in Table 3.

According to the most recent experimental phase diagram of water, upon pressurizing the ambient pressure ice I phase the high-pressure proton ordered phases are expected to occur in the following sequence: ice IX, II, XIII, XIV, XV, and VIII. Salzmann et al. (2011) However, the exact phase boundaries between these phases have not been determined directly from experiment, especially when considering the low temperature regime. Petrenko and Whitworth (2003); Salzmann et al. (2009); sal (); Salzmann et al. (2011) Specifically, between ice I and IX there is no measured phase boundary available and a reasonable choice is to consider the known phase coexistence line between ice I and III, the proton disordered counterpart of ice IX, which appears at 0.1-0.2 GPa. Kell and Whalley (1968); Salzmann et al. (2011) The experimental phase boundaries between ice IX, II, and XIII are also unknown, however, they certainly should appear in the pressure window of 0.2-0.8 GPa since at higher pressures (1.2-1.4 GPa) ice XIV and XV are found to be stabilized. Salzmann et al. (2006, 2009, 2011) The highest density phase, ice VIII, can be found at 1.5 GPa. Salzmann et al. (2011) Now we will discuss how our calculated phase transition pressures compare with the experimental data.

Expt. 0.1–0.2 0.2–0.8 1.2–1.4 1.50
PBE 1.66 3.45 6.08
PBE0 1.32 2.83 4.50
PBE+vdW 0.26 1.25 6.37
PBE0+vdW 0.10 0.89 1.04 4.62
PBE0+vdW 0.26 1.65 3.50
revPBE-vdW 0.78 2.36
optPBE-vdW 0.68 2.22
rPW86-vdW2 0.63 1.25 1.34
Ref. Kell and Whalley, 1968; Ref. Salzmann et al.,2011
Table 3: Comparisons of the calculated and experimental transition pressures. Only positive transition pressures are reported. All pressures are in GPa.

Table 3 shows that the phase transition pressures obtained from PBE are much too high compared to experiment; about an order of magnitude too high for ice IX and 3-4 times too high for ice XIV and VIII. Small improvements arise using PBE0 with 20-30% reductions in the transition pressures. The predicted order in which the ice phases appear (I, IX, XIV, and VIII) with increasing pressure agrees with experiment. However, ice II, XIII, and XV are missing from the PBE (and PBE0) phase diagram, i.e., at no positive external pressure do these phases have the lowest enthalpy (Fig. 4).

With the vdW inclusive functionals the transition pressures are lowered substantially and are in reasonable agreement with experiment (Table 3). For ice IX, XIII, and XIV the transition pressures obtained from PBE0+vdW are within the range of experimental values. However, PBE0+vdW fails to reduce the transition pressure of ice VIII mainly because the relative lattice energy of ice VIII with respect to ice XV (50 meV/HO) does not improve with vdW (Table 1). Inclusion of many-body vdW decreases the energy difference between ice VIII and ice XV by 10 meV/HO and brings the relative change in the equilibrium volume (2.2 Å/HO) into better agreement with experiment (2.4 Å/HO). Both improvements help in reducing the calculated transition pressure of ice VIII to 3.5 GPa which is closer to the experimental pressure of 1.5 GPa. Despite the improvements in the transition pressures both vdW and vdW fail to predict the presence of all of the experimentally characterized ice phases on the phase diagram (Fig. 4).

The phase diagrams obtained with the three vdW-DFs are not particularly impressive either. None of the functionals find the ice I to IX transition at a positive pressure, because they predict ice IX is to be energetically more stable than ice I at zero pressure. However, the predicted transition pressures for the higher density phases (ice XIII and beyond) are in good agreement with experiment, differing by no more than a factor of 2. Interestingly, since ice XIV and XV are isoenergetic with rPW86-vdW2 (Table 1) this functional predicts ice XV to be more stable than ice XIV at all pressures (Fig. 4(f)). For ice VIII all three vdW-DFs reproduce the experimental transition pressure (1.5 GPa) with reasonable accuracy, rPW86-vdW2 being the closest (1.34 GPa) followed by optPBE-vdW (2.22 GPa) and then revPBE-vdW (2.36 GPa).

iii.4 Decomposition of the Exchange and Correlation Contributions to the Lattice Energy

Figure 5: Contributions to the lattice energies of the various ice phases from (a) the exchange energy () [c.f. Eq. 5], (b) the total correlation energy (), and (c) the correlation energy beyond GGA-PBE correlation (). All energies have been calculated on the equilibrium densities obtained from each functional.
revPBE- optPBE- rPW86- PBE+
vdW vdW vdW2 vdW
Isolated HO 58.42 58.33 20.85 43.96
I 61.53 59.30 28.55 48.40
IX 64.31 61.41 33.33 51.36
II 63.88 60.49 33.35 52.24
XIV 64.25 60.60 34.67 53.32
XV 64.67 60.72 35.28 53.86
VIII 64.40 59.19 36.53 55.86
The corresponding experimental value is 45.29 Margoliash and Meath, 1978.
Table 4: The calculated molecular coefficients of water molecules in the various ice phases and an isolated HO molecule. The coefficients are given in HartreeBohr.

GGA, hybrid-GGA, and vdW inclusive functionals lead to varied results for the ice phases considered. In order to shed more light on why this is we have decomposed the contributions from exchange and correlation energies to the lattice energies for all the xc functionals studied. The contribution from the exchange energy () to the lattice energy is obtained by subtracting the exchange energy of isolated HO molecules () from the exchange energy of the ice unit cells () containing HO molecules and can be defined as:


An equivalent definition is used to extract the contribution form the correlation energy () to the lattice energy.

Figs. 5(a) and 5(b) show the variations in and for all ice phases at the equilibrium densities obtained from each xc functional. In general we find that upon going from the low to the high density phases the energetically favorable exchange contribution to the lattice energy decreases, just as the lattice energies do. We also find that the exchange contribution to the lattice energy strongly depends on the equilibrium volumes obtained with the various functionals. Consequently revPBE-vdW predicts the largest volumes and smallest and PBE+vdW the smallest volumes and largest for all phases. Since the hybrid PBE0 exchange yields accurate electrostatic properties (e.g., polarizability, dipole moment, electronic band gap) for the gas and condensed phases of water Hammond et al. (2009); Santra et al. (2009); Wang et al. (2010); Labat et al. (2011) it’s somewhat useful to consider PBE0 as a reference against which we compare other exchange functionals. In this regard PBE and optPBE follow PBE0 fairly closely over the entire pressure range. On the other hand, while rPW86-vdW2 is within 10 meV/HO of the PBE0 value for ice I and IX, it deviates substantially (70 meV/HO more negative) for the higher density ice phases, implying that rPW86 exchange is over stabilized compared to PBE0 exchange for the highest pressure phase. Similarly PBE+vdW and PBE0+vdW are substantially more negative than PBE0 exchange for all phases except ice VIII. However, this is mainly due to the smaller volumes predicted by these approaches compared to PBE0.

The contributions from correlations to the lattice energies, , show why PBE and PBE0 perform so poorly for the high-density phases. Specifically we find that the PBE is nearly constant for all the ice phases, which is in stark contrast to the predictions from the vdW inclusive xc functionals that increases from the low to high density phases (Fig. 5(b)). It is interesting to compare the relative contributions of vdW forces coming from the various vdW inclusive functionals. However, this not straightforward because the correlation energies in vdW and the vdW-DFs contain different terms (c.f. Eqns. 2 and  3). Nonetheless, since in the vdW scheme used here vdW is the correlation energy coming beyond GGA PBE () we have computed a similar quantity from the vdW-DFs by subtracting (with GGA PBE) from . The contribution of this modified non-local correlation energy to the lattice energy of ice is denoted as and is shown in Fig. 5(c). When we examine this term we find that it increases from the low to the high-density ice phases with all vdW inclusive functionals. However, the magnitude of predicted by the different approaches differs significantly. rPW86-vdW2 predicts the smallest , revPBE-vdW and optPBE-vdW the largest, and PBE/PBE0+vdW falls in the middle. Since a major component of is non-local vdW interactions the magnitude of should depend strongly on the vdW coefficients. Indeed, we find that the molecular coefficients calculated on ice (Table 4) with the different functionals correlate well with the relative magnitude of not (b) Previous work showed that the coefficient of an isolated water molecule is 50% too small with rPW86-vdW2 compared to experiment. Vydrov and Voorhis (2010) Here we find the same behavior for the coefficients of water molecules within all ice phases. Compared to all other vdW functionals the molecular coefficients obtained from rPW86-vdW2 are strikingly smaller, 40-50% for ice I and 35-45% for ice VIII.

To sum up, this brief analysis of the exchange and correlation contributions to the lattice energies has revealed that the large reduction in the exchange contribution to the lattice energy upon going from the low- to the high-density phases is compensated for by a growing correlation contribution to the lattice energy from the beyond GGA correlation (). This compensation is obviously found for the vdW inclusive methods but not found for PBE and PBE0 and as a result the lattice energies are underestimated with PBE and PBE0 for the high pressure phases.

Iv Conclusions

We have performed a detailed study on a selection of different ice phases with a range of xc functionals, including some of the recently developed functionals which account for vdW dispersion forces. Whilst we know a lot about the performance of these functionals in the gas phase (in particular on gas phase data sets such as the S22 Jurečka et al. (2006); Burns et al. (2011)) much less is known about how these functionals perform in the condensed phases, which was one of the key motivations for this study. As seen before in the gas phase the vdW inclusive functionals do offer some improvement in performance. This is particularly true for the relative energies of the different phases and as a result the phase transition pressures. However, the functionals tested are far from perfect and none simultaneously yields excellent lattice energies and lattice constants for all phases. Of the schemes considered PBE0+vdW consistently overestimates lattice energies by 50 meV/HO and equilibrium densities by 5%. optPBE-vdW produces densities of ice that are in best agreement (3%) with experiment but the lattice energies are 50 meV/HO too large. revPBE-vdW underestimates densities by 10% and lattice energies by 50 meV/HO. rPW86-vdW2 gives very accurate lattice energies but the densities are underestimated by 8%.

The improved agreement between the experimental and calculated phase transition pressures when using the vdW functionals clearly highlights the importance of accounting for vdW in ice. However, even with vdW inclusive functionals, capturing all of the experimentally characterized ice phases on the water phase diagram is clearly still a major challenge and beyond the capabilities of the methods considered here. Water is well known to provide a stern challenge for DFT, be it water clusters, liquid water and now ice. The fact that several phases of ice are missing from the phase diagram of water is somewhat of a blow to the true predictive ability of the methods considered here, but also a challenge and opportunity for developing and testing new methods.

From this study it is evident that the ice phases considered here are extremely useful in providing a challenging “data set” against which new methods can be tested and proved. It would of course be interesting to see how some of the other vdW inclusive DFT methods developed recently perform on the ice phases. In this respect the already available experimental lattice energies and the matching DMC numbers are valuable references. However, additional DMC data on other phases of ice would certainly be of value as would other vdW inclusive methods e.g., random-phase approximation Lu et al. (2009); Schimka et al. (2010) and second order Møller-Plesset perturbation theory. Grüneis et al. (2010); Hermann and Schwerdtfeger (2008) Finally, we note that the difficulty in predicting ice phases up to only the pressure range 1-2 GPa using GGA, hybrid, and vdW inclusive DFT approaches suggests that caution must be exercised when searching for and predicting new phases of water at yet higher pressures using such functionals. Pickard and Needs (2007); Militzer and Wilson (2010); McMahon (2011); Hermanna et al. (2012); Pickard et al. (2013)

Acknowledgements: B. Santra and R.C. are supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by Department of Energy grant DE-SC0008626. A.M. was supported by the European Research Council (Quantum-CRASS project) and the Royal Society through a Royal Society Wolfson Research Merit Award. Also partly via our membership of the UK’s HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/F067496), this work made use of the computational facilities of HECToR. J.K. is grateful to UCL and the EPSRC for support through the PhD+ scheme.


  • Santra et al. (2007) B. Santra, A. Michaelides, and M. Scheffler, J. Chem. Phys. 127, 184104 (2007).
  • Santra et al. (2008) B. Santra, A. Michaelides, M. Fuchs, A. Tkatchenko, C. Filippi, and M. Scheffler, J. Chem. Phys. 129, 194111 (2008).
  • Santra et al. (2009) B. Santra, A. Michaelides, and M. Scheffler, J. Chem. Phys. 131, 124509 (2009).
  • Dahlke et al. (2008) E. E. Dahlke, R. M. Olson, H. R. Leverentz, and D. G. Truhlar, J. Phys. Chem. A 112, 3976 (2008).
  • Xu and Goddard III (2004) X. Xu and W. A. Goddard III, J. Phys. Chem. A 108, 2305 (2004).
  • Su et al. (2004) J. T. Su, X. Xu, and W. A. Goddard III, J. Phys. Chem. A 108, 10518 (2004).
  • Xantheas (1995) S. S. Xantheas, J. Chem. Phys. 102, 4505 (1995).
  • Kim and Jordan (1994) K. Kim and K. D. Jordan, J. Phys. Chem. 98, 10089 (1994).
  • Anderson and Tschumper (2006) J. A. Anderson and G. S. Tschumper, J. Phys. Chem. A 110, 7268 (2006).
  • Shields and Kirschner (2008) G. C. Shields and K. N. Kirschner, Synthesis and Reactivity in Inorganic, Metal-Organic, and Nano-Metal Chemistry 38, 32 (2008).
  • Csonka et al. (2005) G. I. Csonka, A. Ruzsinszky, and J. P. Perdew, J. Phys. Chem. B 109, 21471 (2005).
  • Dahlke and Truhlar (2005) E. E. Dahlke and D. G. Truhlar, J. Phys. Chem. B 109, 15677 (2005).
  • Gillan et al. (2012) M. J. Gillan, F. R. Manby, M. D. Towler, and D. Alfè, J. Chem. Phys. 136, 244105 (2012).
  • Ireta et al. (2004) J. Ireta, J. Neugebauer, and M. Scheffler, J. Phys. Chem. A 108, 5692 (2004).
  • Tsuzuki and Lüthi (2001) S. Tsuzuki and H. P. Lüthi, J. Chem. Phys. 114, 3949 (2001).
  • Novoa and Sosa (1995) J. J. Novoa and C. Sosa, J. Phys. Chem. 99, 15837 (1995).
  • Hammond et al. (2009) J. R. Hammond, N. Govind, K. Kowalski, J. Autschbach, and S. S. Xantheas, J. Chem. Phys. 131, 214103 (2009).
  • Wang et al. (2010) F.-F. Wang, G. Jenness, W. A. Al-Saidi, and K. D. Jordan, J. Chem. Phys. 132, 134303 (2010).
  • Grossman et al. (2004) J. C. Grossman, E. Schwegler, E. W. Draeger, F. Gygi, and G. Galli, J. Chem. Phys. 120, 300 (2004).
  • Schwegler et al. (2004) E. Schwegler, J. C. Grossman, F. Gygi, and G. Galli, J. Chem. Phys. 121, 5400 (2004).
  • Fernández-Serra and Artacho (2004) M. V. Fernández-Serra and E. Artacho, J. Chem. Phys. 121, 11136 (2004).
  • McGrath et al. (2005) M. J. McGrath, J. I. Siepmann, I.-F. W. Kuo, C. J. Mundy, J. VandeVondele, J. Hutter, F. Mohamed, and M. Krack, ChemPhysChem 6, 1894 (2005).
  • McGrath et al. (2006a) M. J. McGrath, J. I. Siepmann, I.-F. W. Kuo, and C. Mundy, Molec. Phys. 104, 3619 (2006a).
  • McGrath et al. (2006b) M. J. McGrath, J. I. Siepmann, I.-F. W. Kuo, C. J. Mundy, J. VandeVondele, J. Hutter, F. Mohamed, and M. Krack, J. Phys. Chem. A 110, 640 (2006b).
  • Lee and Tuckerman (2006) H.-S. Lee and M. E. Tuckerman, J. Chem. Phys. 125, 154507 (2006).
  • Lee and Tuckerman (2007) H.-S. Lee and M. E. Tuckerman, J. Chem. Phys. 126, 164501 (2007).
  • Todorova et al. (2006) T. Todorova, A. P. Seitsonen, J. Hutter, I.-F. W. Kuo, and C. J. Mundy, J. Phys. Chem. B 110, 3685 (2006).
  • VandeVondele et al. (2005) J. VandeVondele, F. Mohamed, M. Krack, J. Hutter, M. Sprik, and M. Parrinello, J. Chem. Phys. 122, 014515 (2005).
  • Guidon et al. (2008) M. Guidon, F. Schiffmann, J. Hutter, and J. VandeVondele, J. Chem. Phys. 128, 214104 (2008).
  • Kühne et al. (2009) T. D. Kühne, M. Krack, and M. Parrinello, J. Chem. Theory Comput. 5, 235 (2009).
  • Sit and Marzari (2005) P. H.-L. Sit and N. Marzari, J. Chem. Phys. 122, 204510 (2005).
  • Fernández-Serra et al. (2005) M. V. Fernández-Serra, G. Ferlat, and E. Artacho, Molecular Simualtion 31, 361 (2005).
  • Asthagiri et al. (2003) D. Asthagiri, L. R. Pratt, and J. D. Kress, Phys. Rev. E 68, 041505 (2003).
  • Mattson and Mattson (2009) A. E. Mattson and T. R. Mattson, J. Chem. Theory Comput. 5, 887 (2009).
  • Kuo et al. (2004) I.-F. W. Kuo, C. J. Mundy, M. J. McGrath, J. I. Siepmann, J. VandeVondele, M. Sprik, J. Hutter, B. Chen, M. L. Klein, F. Mohamed, et al., J. Phys. Chem. B 108, 12990 (2004).
  • Sprik et al. (1996) M. Sprik, J. Hutter, and M. Parrinello, J. Chem. Phys. 105, 1142 (1996).
  • Silvestrelli and Parrinello (1999) P. L. Silvestrelli and M. Parrinello, J. Chem. Phys. 111, 3572 (1999).
  • Yoo et al. (2009) S. Yoo, X. C. Zeng, and S. S. Xantheas, J. Chem. Phys. 130, 221102 (2009).
  • Zhang et al. (2011a) C. Zhang, D. Donadio, F. Gygi, and G. Galli, J. Chem. Theory Comput. 7, 1443 (2011a).
  • Chen et al. (2003) B. Chen, I. Ivanov, M. L. Klein, and M. Parrinello, Phys. Rev. Lett. 91, 215503 (2003).
  • Morrone and Car (2008) J. A. Morrone and R. Car, Phys. Rev. Lett. 101, 017801 (2008).
  • Laasonen et al. (1993) K. Laasonen, M. Sprik, M. Parrinello, and R. Car, J. Chem. Phys. 99, 9080 (1993).
  • Akin-Ojo and Wang (2011) O. Akin-Ojo and F. Wang, Chem. Phys. Lett. 513, 59 (2011).
  • Alfè et al. (2013) D. Alfè, A. P. Bartók, G. Csányi, and M. J. Gillan, J. Chem. Phys. 138, 221102 (2013).
  • Santra et al. (2011) B. Santra, J. Klimeš, D. Alfè, A. Tkatchenko, B. Slater, A. Michaelides, R. Car, and M. Scheffler, Phys. Rev. Lett. 107, 185701 (2011).
  • Murray and Galli (2012) E. D. Murray and G. Galli, Phys. Rev. Lett. 108, 105502 (2012).
  • Hamada (2010) I. Hamada, J. Chem. Phys. 133, 214503 (2010).
  • Kolb and Thonhauser (2011) B. Kolb and T. Thonhauser, Phys. Rev. B 84, 045116 (2011).
  • Kambara et al. (2012) O. Kambara, K. Takahashi, M. Hayashi, and J.-L. Kuo, Phys. Chem. Chem. Phys. 14, 11484 (2012).
  • Fang et al. (2013) Y. Fang, B. Xiao, J. Tao, J. Sun, and J. P. Perdew, Phys. Rev. B 87, 214101 (2013).
  • Gillan et al. (2013) M. J. Gillan, D. Alfè, A. P. Bartók, and G. Csányi, arXiv:1303.0751 (2013), http://arxiv.org/abs/1303.0751.
  • Pamuk et al. (2012) B. Pamuk, J. M. Soler, R. Ramírez, C. P. Herrero, P. Stephens, P. B. Allen, and M.-V. Fernández-Serra, Phys. Rev. Lett. 108, 193003 (2012).
  • Lee et al. (1992) C. Lee, D. Vanderbilt, K. Laasonen, R. Car, and M. Parrinello, Phys. Rev. Lett. 69, 462 (1992).
  • Hamann (1997) D. R. Hamann, Phys. Rev. B 55, R10157 (1997).
  • Singer et al. (2005) S. J. Singer, J.-L. Kuo, T. K. Hirsch, C. Knight, L. Ojamäe, and M. L. Klein, Phys. Rev. Lett. 94, 135701 (2005).
  • Tribello et al. (2006) G. A. Tribello, B. Slater, and C. G. Salzmann, J. Am. Chem. Soc. 128, 12594 (2006).
  • de Koning et al. (2006) M. de Koning, A. Antonelli, A. J. R. da Silva, and A. Fazzio, Phys. Rev. Lett. 96, 075501 (2006).
  • Feibelman (2008) P. J. Feibelman, Phys. Chem. Chem. Phys. 10, 4688 (2008).
  • Hermann and Schwerdtfeger (2008) A. Hermann and P. Schwerdtfeger, Phys. Rev. Lett. 101, 183005 (2008).
  • Pan et al. (2008) D. Pan, L.-M. Liu, G. A. Tribello, B. Slater, A. Michaelides, and E. Wang, Phys. Rev. Lett. 101, 155703 (2008).
  • Erba et al. (2009) A. Erba, S. Casassa, L. Maschio, and C. Pisani, J. Phys. Chem. B 113, 2347 (2009).
  • Morrone et al. (2009) J. A. Morrone, L. Lin, and R. Car, J. Chem. Phys. 130, 204511 (2009).
  • Militzer and Wilson (2010) B. Militzer and H. F. Wilson, Phys. Rev. Lett. 105, 195701 (2010).
  • Labat et al. (2011) F. Labat, C. Pouchan, C. Adamo, and G. E. Scuseria, J. Comput. Chem. 32, 2177 (2011).
  • Lin et al. (2009) I.-C. Lin, A. P. Seitsonen, M. D. Coutinho-Neto, I. Tavernelli, and U. Rothlisberger, J. Phys. Chem. B 113, 1127 (2009).
  • Schmidt et al. (2009) J. Schmidt, J. VandeVondele, I.-F. W. Kuo, D. Sebastiani, J. I. Siepmann, J. Hutter, and C. J. Mundy, J. Phys. Chem. B 113, 11959 (2009).
  • Wang et al. (2011) J. Wang, G. Román-Pérez, J. M. Soler, E. Artacho, and M.-V. Fernández-Serra, J. Chem. Phys. 134, 024516 (2011).
  • Møgelhøj et al. (2011) A. Møgelhøj, A. K. Kelkkanen, K. T. Wikfeldt, J. Schiøtz, J. J. Mortensen, L. G. M. Pettersson, B. I. Lundqvist, K. W. Jacobsen, A. Nilsson, and J. K. Nørskov, J. Phys. Chem. B 115, 14149 (2011).
  • Zhang et al. (2011b) C. Zhang, J. Wu, G. Galli, and F. Gygi, J. Chem. Theory Comput. 7, 3054 (2011b).
  • Jonchiere et al. (2011) R. Jonchiere, A. P. Seitsonen, G. Ferlat, A. M. Saitta, and R. Vuilleumier, J. Chem. Phys. 135, 154503 (2011).
  • Yoo and Xantheas (2012) S. Yoo and S. S. Xantheas, J. Chem. Phys. 134, 121105 (2012).
  • Ma et al. (2012) Z. Ma, Y. Zhang, and M. E. Tuckerman, J. Chem. Phys. 137, 044506 (2012).
  • (73) Z. Li, Improving Ab Initio Molecular Dynamics of Liquid Water, Ph.D. Thesis, Princeton University, 2012, http://www.princeton.edu/physics/graduate-program/theses/Thes%is_Zhaofeng-Li.pdf.
  • Dion et al. (2004) M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
  • Tkatchenko and Scheffler (2009) A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
  • Klimeš et al. (2010) J. Klimeš, D. R. Bowler, and A. Michaelides, J. Phys.: Cond. Matt. 22, 022201 (2010).
  • von Lilienfeld et al. (2004) O. A. von Lilienfeld, I. Tavernelli, U. Rothlisberger, and D. Sebastiani, Phys. Rev. Lett. 93, 153004 (2004).
  • Silvestrelli (2008) P. L. Silvestrelli, Phys. Rev. Lett. 100, 053002 (2008).
  • Sato and Nakai (2009) T. Sato and H. Nakai, J. Chem. Phys. 131, 224104 (2009).
  • Grimme et al. (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  • Vydrov and Voorhis (2009) O. A. Vydrov and T. V. Voorhis, Phys. Rev. Lett. 103, 063004 (2009).
  • Becke and Johnson (2005) A. D. Becke and E. R. Johnson, J. Chem. Phys. 123, 154101 (2005).
  • Tao et al. (2012) J. Tao, J. P. Perdew, and A. Ruzsinszky, Proc. Natl. Acad. Sci. (USA) 109, 18 (2012).
  • Petrenko and Whitworth (2003) V. Petrenko and R. W. Whitworth, Physics of Ice (Oxford University Press, New York, 2003).
  • Salzmann et al. (2009) C. G. Salzmann, P. G. Radaelli, E. Mayer, and J. L. Finney, Phys. Rev. Lett. 103, 105701 (2009).
  • (86) C. G. Salzmann, P. G. Radaelli, A. Hallbrucker, E. Mayer, and J. L. Finney, in Physics and Chemistry of Ice, ed. W. F. Kuhs, The Royal Society of Chemistry, Cambridge, 2007, pp. 521–528.
  • Salzmann et al. (2011) C. G. Salzmann, P. G. Radaelli, B. Slater, and J. L. Finney, Phys. Chem. Chem. Phys. 13, 18468 (2011).
  • Tkatchenko et al. (2012) A. Tkatchenko, R. A. DiStasio Jr., R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012).
  • Lee et al. (2010) K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, Phys. Rev. B 82, 081101(R) (2010).
  • Langreth et al. (2009) D. C. Langreth et al., J. Phys.: Condens. Matter 21, 084203 (2009).
  • Tkatchenko et al. (2010) A. Tkatchenko, L. Romaner, O. T. Hofmann, E. Zojer, and C. Ambrosch-Draxl, MRS Bull. 35, 435 (2010).
  • Klimeš et al. (2011) J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011).
  • DiStasio Jr. et al. (2012) R. A. DiStasio Jr., O. A. von Lilienfeld, and A. Tkatchenko, Proc. Natl. Acad. Sci. (USA) 109, 14791 (2012).
  • Klimeš and Michaelides (2012) J. Klimeš and A. Michaelides, J. Chem. Phys. 137, 120901 (2012).
  • Reilly and Tkatchenko (2013) A. M. Reilly and A. Tkatchenko, J. Phys. Chem. Lett. 4, 1028 (2013).
  • Kelkkanen et al. (2009) A. K. Kelkkanen, B. I. Lundqvist, and J. K. Nørskov, J. Chem. Phys. 131, 046102 (2009).
  • Line and Whitworth (1996) C. M. B. Line and R. W. Whitworth, J. Chem. Phys. 104, 10008 (1996).
  • Londono et al. (1993) J. D. Londono, W. F. Kuhs, and J. L. Finney, J. Chem. Phys. 98, 4878 (1993).
  • Lobban et al. (2002) C. Lobban, J. L. Finney, and W. F. Kuhs, J. Chem. Phys. 117, 3928 (2002).
  • Salzmann et al. (2006) C. G. Salzmann, P. G. Radaelli, A. Hallbrucker, E. Mayer, and J. L. Finney, Science 311, 1758 (2006).
  • Kuhs et al. (1984) W. F. Kuhs, J. L. Finney, C. Vettier, and D. V. Bliss, J. Chem. Phys. 81, 3612 (1984).
  • Pan et al. (2010) D. Pan, L.-M. Liu, G. A. Tribello, B. Slater, A. Michaelides, and E. Wang, J. Phys.: Condensed Matter 22 22, 074203 (2010).
  • Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  • (104) F. D. Murnaghan, Proc. Natl. Acad. Sci. USA 30, 244 (1944).
  • (105) B. Santra, Density-Functional Theory Exchange-Correlation Functionals for Hydrogen Bonds in Water, Ph.D. Thesis, Fritz-Haber-Institut der Max-Planck Gesellschaft, TU-Berlin, 2010, http://opus4.kobv.de/opus4-tuberlin/frontdoor/index/index/doc%Id/2723.
  • Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  • Zhang and Yang (1998) Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998).
  • Jurečka et al. (2006) P. Jurečka, J. Šponer, J. Černý, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006).
  • Murray et al. (2009) E. D. Murray, K. Lee, and D. C. Langreth, J. Chem. Theory Comput. 5, 2754 (2009).
  • Blum et al. (2009) V. Blum et al., Comp. Phys. Comm. 180, 2175 (2009).
  • Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  • Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  • Román-Pérez and Soler (2009) G. Román-Pérez and J. M. Soler, Phys. Rev. Lett. 103, 096102 (2009).
  • Whalley (1984) E. Whalley, J. Chem. Phys. 81, 4087 (1984).
  • Röttger et al. (1994) K. Röttger, A. Endriss, J. Ihringer, S. Doyle, and W. F. Kuhs, Acta Cryst. B50, 644 (1994).
  • Carrasco et al. (2011) J. Carrasco, B. Santra, J. Klimeš, and A. Michaelides, Phys. Rev. Lett. 106, 026101 (2011).
  • Gulans et al. (2009) A. Gulans, M. J. Puska, and R. M. Nieminen, Phys. Rev. B 79, 201105(R) (2009).
  • not (a) Frequency calculations have not been performed for any of the PBE0 based functionals. For PBE0 we have estimated the ZPE corrected volumes by considering the volume changes found with PBE. For PBE0+vdW and PBE0+vdW ZPE corrected volumes are estimated with the volume changes obtained with PBE+vdW.
  • Kell and Whalley (1968) G. S. Kell and E. Whalley, J. Chem. Phys. 48, 2359 (1968).
  • not (b) In case of vdW molecular coefficients are calculated by using Eq. 10 in Ref. Tkatchenko and Scheffler (2009) and for the vdW-DFs we have used the long distance behavior of the kernel as given by Eq. 17 in Ref. Dion et al. (2004) and in Eq. 23.5 in Ref. Vydrov and Van Voorhis (2012). The vdW-DF coefficients vary only within 1% when using different PAW potentials or norm-conserving pseudopotentials, however, they depend significantly on the functional employed to calculate the electron density. For example, the coefficient of an isolated water molecule is reduced to 45.5 a.u. when Hartree-Fock electron density is used. Likewise, the vdW-DF coefficients of an isolated water molecule obtained here are 20-24% larger than the values reported in Ref. Vydrov and Voorhis (2010) where electron densities were calculated using a so-called long-range corrected hybrid functional which is expected to give densities similar to those given by Hartree-Fock.
  • Vydrov and Voorhis (2010) O. A. Vydrov and T. V. Voorhis, Phys. Rev. A 81, 062708 (2010).
  • Burns et al. (2011) L. A. Burns, Álvaro Vázquez-Mayagoitia, B. G. Sumpter, and C. D. Sherrill, J. Chem. Phys. 134, 084107 (2011).
  • Lu et al. (2009) D. Lu, Y. Li, D. Rocca, and G. Galli, Phys. Rev. Lett. 102, 206411 (2009).
  • Schimka et al. (2010) L. Schimka, J. Harl, A. Stroppa, A. Grüneis, M. Marsman, F. Mittendorfer, and G. Kresse, Nature Mater. 9, 741 (2010).
  • Grüneis et al. (2010) A. Grüneis, M. Marsman, and G. Kresse, J. Chem. Phys. 133, 074107 (2010).
  • Pickard and Needs (2007) C. J. Pickard and R. J. Needs, J. Chem. Phys. 127, 244503 (2007).
  • McMahon (2011) J. M. McMahon, Phys. Rev. B 84, 220104(R) (2011).
  • Hermanna et al. (2012) A. Hermanna, N. W. Ashcroft, and R. Hoffmann, Proc. Natl. Acad. Sci. (USA) 109, 745 (2012).
  • Pickard et al. (2013) C. J. Pickard, M. Martinez-Canales, and R. J. Needs, Phys. Rev. Lett. 110, 245701 (2013).
  • Vydrov and Van Voorhis (2012) O. A. Vydrov and T. Van Voorhis, in Fundamentals of Time-Dependent Density Functional Theory, edited by M. A. L. Marques, N. Maitra, F. Nogueira, E. K. U. Gross, and A. Rubio (Springer, Berlin, 2012).
  • (131) Supplementary material provides coordinates of all ice phases optimized with all xc functionals considered here. It can be obtained via https://docs.google.com/file/d/0B_RP_vP8_oA5SzItT0ZSRW1Pa1U/edit.
  • Margoliash and Meath (1978) D. J. Margoliash and W. J. Meath, J. Chem. Phys. 68, 1426 (1978).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description