Benchmarking the performance of Density Functional Theory and Point Charge Force Fields in their Description of sI Methane Hydrate against Diffusion Monte Carlo

High quality reference data from diffusion Monte Carlo calculations are presented for bulk sI methane hydrate, a complex crystal exhibiting both hydrogen-bond and dispersion dominated interactions. The performance of some commonly used exchange-correlation functionals and all-atom point charge force fields is evaluated. Our results show that none of the exchange-correlation functionals tested are sufficient to describe both the energetics and the structure of methane hydrate accurately, whilst the point charge force fields perform badly in their description of the cohesive energy but fair well for the dissociation energetics. By comparing to ice Ih, we show that a good prediction of the volume and cohesive energies for the hydrate relies primarily on an accurate description of the hydrogen bonded water framework, but that to correctly predict stability of the hydrate with respect to dissociation to ice Ih and methane gas, accuracy in the water-methane interaction is also required. Our results highlight the difficulty that density functional theory faces in describing both the hydrogen bonded water framework and the dispersion bound methane.


I. INTRODUCTION
The clathrate hydrates of natural gases -crystalline compounds in which gas is dissolved in a host framework of water molecules -are important to a wide variety of applications across the energy and climate sciences. For example, the fact that one volume of hydrate can generate up to 180 volumes of gas upon dissociation at standard temperature and pressure, whilst only 15% of the recovered energy is required for dissociation, means that hydrate reservoirs are a potential untapped energy resource. 1 Even though there remains uncertainty in the total amount of hydrated gas on Earth, there is a consensus that this amount exceeds conventional gas reserves by at least an order of magnitude. 2 Perhaps a more pressing issue is that hydrates also pose a severe problem for flow assurance in oil and gas pipelines: if the mixed phases of water and natural gas are allowed to cool, hydrates may form and block the line, causing production to stall. As readily available oil and gas reserves become depleted, and the need for extraction from deeper reservoirs increases, the consequences of hydrate formation are becoming more severe. Although chemicals for inhibiting hydrate formation exist, they have generally been found on a trial-and-error basis, with little understanding of how they work at the molecular scale. This state of affairs has arisen from the fact that we have little a) Electronic mail: angelos.michaelides@ucl.ac.uk knowledge of the fundamental mechanisms that underlie hydrate formation. Consequently, computer simulation has been used in recent years in attempts to improve our molecular level understanding of hydrate formation. [3][4][5][6][7] It is important, therefore, to understand both the molecular interactions present in condensed phase gas hydrates, and the performance of current approximations used to describe these interactions.
By far the most commonly used electronic structure method for investigating condensed phase systems is density functional theory (DFT) (a recent review of DFT and the current challenges it faces is given in Ref. 8). Despite incredible success in its application to a wide variety of systems, DFT has a number of limitations. Of particular relevance to gas hydrates is the known deficiency of the local density approximation (LDA) and generalized gradient approximation (GGA) varieties of exchange-correlation (xc) functionals to account for van der Waals (vdW) dispersion interactions. Incorporating an accurate description of vdW interactions into density functional theory is a very active research area, with recent developments including Grimme's dispersion correction methods, 9,10 the Tkatchenko-Scheffler scheme 11 and the fully self-consistent vdW-DF method of Dion et al. 12 and its various derivatives. 13,14 For a recent overview of these and other methods to incorporate vdW interactions into DFT see Ref. 15. Understanding the contribution of vdW to the bonding in solids is an important issue and there is a need to understand the strengths and weaknesses of various vdW-inclusive methods in order to improve the performance of DFT. Recent work has shown that vdW-inclusive DFT methods offer a systematic improvement over GGA functionals in describing the phase behaviour of ice. 16,17 Like ice, gas hydrates also have an extended hydrogen bonded network of water molecules, but unlike ice, they contain cavities that gas molecules can occupy. Natural gas hydrates therefore offer the opportunity to test the ability of vdW-inclusive methods to simultaneously describe both the hydrogen bonded water network and the predominantly dispersion bound watergas interaction. As well as DFT, force fields (FFs) are often used to investigate gas hydrates, [3][4][5]18,19 especially when long time and length scales are required, such as in the study of nucleation processes.
Evaluating the performance of techniques such as DFT or FFs requires high-quality reference data to compare to -something that is lacking for gas hydrates in the condensed phase. For example, previous DFT studies 20,21 have evaluated the performance of the chosen xc functional through comparison to experiment or quantum chemical methods on isolated clusters. However, various issues can arise when validating the performance of DFT to experiment, such as temperature/pressure, non-stoichiometry and quantum nuclear effects. Furthermore, although a source of valuable information for understanding the nature of interactions, comparison to isolated clusters (to which accurate quantum chemical methods are generally limited to) does not directly tell us how DFT methods are performing for the condensed phase. The tendency to validate FFs used in molecular dynamics or Monte Carlo simulation against experiment is even greater than it is for DFT. One method that has been shown to provide accurate energies for condensed phase water systems is diffusion Monte Carlo (DMC). DMC can be applied to a range of systems, both isolated and periodic, 22,23 has mild scaling behaviour 24,25 and has rapid and automatic basis-set convergence. 26 DMC has also been shown to favour well in comparison to CCSD(T) -the so-called 'gold-standard' quantum chemical method -for calculations on the water dimer and other small water clusters. [27][28][29][30] It also gives a good description of the relative energies of different ice phases 16,17 and has recently been shown to achieve subchemical accuracy for non-covalent interactions in the gas phase. 31 We therefore have confidence that DMC can be used to obtain accurate reference data for periodic gas hydrate crystals. Specifically, we will use DMC to calculate accurate data for the energetics of a methane hydrate crystal.
In this study, we compare the performance of a number of different xc functionals, ranging from the LDA and PBE 32 levels of approximation, in which vdW interactions are not accounted for, to a variety of dispersioncorrected functionals, namely: an empirical correction scheme from Grimme (PBE-D2); 9 the method developed by Tkatchenko and Scheffler (PBE-vdW TS ), 11 which like PBE-D2 involves an explicit summation of pairwise vdW dispersion interactions over all atom pairs, but differs in that the vdW C 6 coefficients are themselves functionals of the electron density; and a number of functionals from the vdW-DF family. In particular, we consider the original vdW-DF of Dion et al. and the modified versions of Klimeš et al., 13,14 in which the exchange functional is changed from that of revPBE, to 'optPBE', 'optB88' and 'optB86b'. These modified versions of vdW-DF have been shown to offer good performance for a wide range of systems. 13,14,33,34 Throughout the rest of the paper, the original vdW-DF of Dion et al. will be referred to as 'revPBE-vdW' with the term 'vdW-DF' used when referring to the class of functionals. We also present results using the OPLS-AA 35 potential for methane and the TIP4P-2005 36 and TIP4P-ICE 37 potentials for water. Details of these FFs are given in the Supporting Information (SI), 38 but key features of these potentials are that they are all atomic, point charge and have Lennard-Jones sites located on the carbon and oxygen atoms. The TIP4P-2005 and TIP4P-ICE potentials are rigid, whereas the OPLS-AA potential is flexible. We have also investigated the two water potentials in combination with a number of different methane potentials, but as these yield similar results to OPLS-AA, they have been omitted from the main article for clarity and are included in the SI. Although this is clearly not an exhaustive list of possible xc functionals and FFs available, our test set nevertheless is adequate to highlight the main strengths and weaknesses of these types of methods in describing hydrogen-bond plus dispersion bound systems such as methane hydrate.
In the following sections, we will compare the results of the above mentioned xc correlation functionals and force fields to DMC in their prediction of the bulk properties of sI methane hydrate. We will specifically look at the cohesive energy of the hydrate crystal, the binding energy of the methane to the water framework and the dissociation energy of the hydrate crystal to ice I h and methane vapour. We will see that none of the methods give a particularly satisfactory description of bulk sI methane hydrate and that in instances of apparent agreement, this is due to a fortuitous cancellation of errors.

II. COMPUTATIONAL SETUP
DFT calculations were performed using VASP 5.3.2, 39-41 a periodic plane-wave basis set code. 42 Calculations with the vdW-DFs have been carried out self-consistently using the scheme of Román-Pérez and Soler, 43 as implemented in VASP by Klimeš et al.. 14 Projector-augmentedwaves (PAW) potentials 44 have been used, with LDAbased potentials used for the LDA calculations and PBEbased potentials used for all other calculations. All results reported here used the 'standard' PAW potentials supplied with VASP and a plane-wave cutoff of 600 eV (these PAWs have been optimised for a plane-wave basis cutoff ≥ 400 eV). A Γ-centred 2 × 2 × 2 Monkhorst-Pack k-point mesh 45 per unit cell was used for calculations of bulk sI methane hydrate, whereas calculations con-cerning isolated molecules were performed at the Γ-point only, in a cubic simulation cell of volume 20 × 20 × 20Å 3 . The structures for bulk sI methane hydrate were taken from the work of Lenz and Ojmäe 19 and optimised using the conjugate gradient geometry optimiser until forces on all atoms were below 0.02 eV/Å. Wave functions were converged to within 1 × 10 −8 eV. For calculations concerning ice I h , we have made use of the same protonordered 12 molecule ice I h unit cell structures as those in Ref. 17, with a Γ-centred 2 × 2 × 2 Monkhorst-Pack k-point mesh per unit cell used. All other settings were identical to those used for the hydrate calculations.
All quantum Monte Carlo calculations were performed using version 2.12.1 of the Cambridge CASINO code. 25 DMC simulations for 178-atom simulation cells were performed using conventional Slater-Jastrow trial wave functions with a Jastrow factor containing electron-nucleus, electron-electron, and electron-nucleus electron terms, 46 each of which depends on variational parameters determined by a combination of variance-and energyminimization. The orbitals in the determinantal part of the trial wave function were generated from DFT calculations performed by the PWSCF component of the Quantum Espresso package; 47 these Γ-point DFT calculations were done using the PBE xc functional and a 300 Ry (4082 eV) plane-wave cutoff. The same structures from Lenz and Ojamäe used for the VASP calculations were first optimised with these settings in the PWSCF component of the Quantum Espresso package. As is standard practice, the plane-wave orbitals were re-expressed in B-splines 26 for the DMC simulations. Dirac-Fock pseudopotentials specifically developed for use in QMC were used. 48,49 Coulomb finite size effects were accounted for using the 'structure factor' method described in Refs. 50 and 51 (though we could equally well have used the Modified Periodic Coulomb (MPC) interaction defined in Refs 52 and 53, which we checked gave essentially the same results).
FF calculations were performed using the GROMACS 4.5.5 simulation package. 54 Long range electrostatics were treated with the particle-mesh Ewald method 55,56 with a grid spacing of 1Å used for the fast Fourier transform (fourth order interpolation was also used) and a real space cut-off of 9Å. Lennard-Jones interactions were truncated after 9Å with tail corrections applied. The calculations were also performed without the tail corrections and results from these have been included in the SI (any effect of the tail corrections does not alter the conclusions presented in the main paper). The L-BFGS algorithm 57,58 was used to optimise the geometries, with the SETTLES algorithm 59 used to constrain the water geometry. All geometries were converged to within 1.05 × 10 −6 eV/Å.

III. RESULTS AND DISCUSSION
Gas hydrates come in three main crystal forms -structures I, II and H (sI, sII and sH, respectively). Methane hydrate is generally found in the sI form, although the sII and sH forms have been reported under very high pressure (above 250 MPa and ca. 1 GPa respectively). 60,61 In the sI hydrate, the water molecules form a hydrogen bonded network that gives rise to two types of cavities: a twelve-sided pentagonal dodecahedron (often denoted as 5 12 ); and a 14-sided tetrakaidecahedron (denoted as 5 12 6 2 , owing to the fact that it consists of 12 pentagonal and 2 hexagonal faces). In stoichiometric sI hydrate, the methane molecules singly occupy each cavity. The cubic unit cell consists of two 5 12 and six 5 12 6 2 cages and has the chemical formula 46H 2 O · 8CH 4 . The sI methane hydrate structure is shown in Fig. 1(b) and a comprehensive overview of the sI, sII and sH hydrate structures can be found in Ref. 1.
We have computed the cohesive energy per water molecule: of the bulk sI methane hydrate unit cell for a variety of unit cell volumes, maintaining a cubic simulation cell. In Eq. 1, E sI (a) is the total energy of bulk sI methane hydrate with lattice constant a, whilst E H2O and E CH 4 are energies of the isolated water and methane molecules, respectively. We did this first using DMC. By fitting ∆E sI coh (a) to Murnaghan's equation of state, 62,63 we are able to determine the equilibrium lattice constant a 0 and cohesive energy ∆E sI coh (a 0 ). These results are presented in Fig. 1(a), where we can see that the equilibrium lattice constant is estimated to be 11.83 ± 0.02Å and the cohesive energy is −632 ± 1 meV/H 2 O. The DMC lattice constant compares well to the low temperature neutron scattering data of Davidson et al. 64 (11.77 ± 0.01Å, CH 4 /D 2 O at 5.2 K) and Gutt et al. 65 (11.821 ± 0.001Å, CD 4 /D 2 O at 2 K).
We have also computed the variation of the cohesive energy with lattice constant for each of the DFT xc functionals and FFs discussed in Section I. In these calculations, all atoms were allowed to relax independently (with the constraint of rigid water molecules for the FFs). The results of these calculations are presented in Figure 2 and Table I (the LDA results have been excluded from Fig. 2 for clarity). Although all of the examined DFT xc functionals overbind the hydrate crystal, there is considerable variety amongst the DFT results, with optPBE-, optB88-, optB86b-vdW, PBE-D2 and PBE-vdW TS significantly overbinding the hydrate crystal, whilst PBE and revPBE-vdW yield cohesive energies in better agreement with DMC. Despite having the best agreement with DMC for the cohesive energy (within 1%), revPBE-vdW does, however, predict a lattice constant that is 1.9-2.3% too large. On the other hand, although it significantly overbinds the crystal, optPBE-vdW predicts a structure in decent agreement with DMC (0.5-0.9% too small), whereas optB88-and optB86b-vdW yield lattice constants that are too short by 2.3-2.6% and 2.5-2.9%, respectively. PBE-D2 and PBE-vdW TS also strongly overbind the hydrate crystal and predict lattice constants that are too small by 3.0% or worse. What is perhaps surprising is that PBE, which fails to account for vdW interactions entirely, is yielding reasonable results not only for the structure (0.6-0.9% smaller than DMC), but also for the energetics. In fact, not only does PBE predict an equilibrium cohesive energy in reasonable agreement with the DMC result, it actually slightly overbinds the crystal by 2.1-2.4%. The force fields, OPLS-AA/TIP4P-2005 and OPLS-AA/TIP4P-ICE, overbind by 8.2-8.6% and 18.0-18.4 % respectively, although their predicted structures are in decent agreement with the reference data, with their predicted lattice constants differing from DMC by less than 1.1%.
From the results for ∆E sI coh (a) presented in Fig. 2 it would be tempting to conclude that PBE gives a satisfactory description of bulk sI methane hydrate. Given the well known problem that GGA functionals do not account for dispersion interactions, however, the fact that PBE slightly overbinds the hydrate seems almost paradoxical. Furthermore, the overly repulsive nature of revPBE exchange at short separations has been shown to lead to lattice constants that are too long and cohesive energies that are too weak in hydrogen bonded systems such as ice. 17 Indeed, we do obtain a lattice constant with revPBE-vdW that is 1.7-2.2% too large, but why then, is the cohesive energy for sI methane hydrate slightly too strong with this functional? To better understand these results, we have decomposed the total cohesive energy into contributions arising from the methane binding to and the cohesive energy of the empty hydrate: where E empty (a 0 ) is the energy of the hydrate unit cell with no methane present, calculated without further relaxation of the water molecules (i.e. the water molecules are 'frozen' in the position they assume in the bulk hydrate). For DMC, both ∆E CH 4 and ∆E empty coh have been calculated at the experimental lattice constant 64 a = 11.77Å. These results are presented in Fig. 3 and Table I, with DMC providing reference values of ∆E CH 4 = −241±15 meV/CH 4 and ∆E empty The origin of PBE's seemingly good description of bulk sI hydrate now becomes apparent: the lack of vdW interactions means there is no binding between the methane and the water (in fact ∆E CH 4 is slightly positive), but this is compensated for by an overbinding of the hydrogen bonded water framework. Although the overbinding of the water framework is small on a per molecule basis, water and methane exist in a ratio of 23:4 in the stoichiometric hydrate, meaning that small errors in describing the water-water interactions are much amplified compared to the apparently larger errors in the methane binding energy. From Fig. 3 we can also see that LDA's severe overbinding occurs principally from its description of the water framework (∆E empty,LDA coh − ∆E empty,DMC coh = −531 meV/H 2 O), although it is worth noting that it also overbinds the methane to the water framework by 85 meV/CH 4 . The ability of LDA to bind van der Waals systems (such as CH 4 in a H 2 O cage) has been observed before; 66,67 this is known to be fortuitous because, by its nature, LDA relies on a local description of exchange and correlation and does not account for non-local interactions. Turning our attention to the dispersion-corrected functionals it is clear that, with the exception of PBE-D2, they all over-correct the neglect of vdW interactions by the GGA functional, yielding methane binding energies that are too strong by 138-262 meV/CH 4 . It is also clear that the better agreement of the cohesive energy obtained with revPBE-vdW compared to the other dispersion-corrected functionals is due to an underbinding of the water framework (consistent with results obtained for bulk ice I h 17 ) that offsets a strong overbinding of the methane. In the case of the other dispersioncorrected functionals, as well as predicting methane binding energies that are too exothermic, they also overbind the water framework by 83-130 meV/H 2 O. The source of overbinding for the FFs occurs almost exclusively in the water framework, with both FFs presented here yielding good agreement for ∆E CH 4 .   . Apart from PBE-D2, all dispersion-corrected density functionals severely overbind methane to the empty hydrate. Similarly, all density functionals overbind the water framework, with the exception of revPBE-vdW. PBE, which does not account for vdW interactions, fails to predict methane binding to the empty hydrate structure. The force fields yield good values for ∆ECH 4 , but like the DFT methods, they overbind the water framework.
Due to the high water content of sI methane hydrate, it is convenient to compare the performance of the xc functionals and FFs for the hydrate to ice I h . We choose ice I h rather than any of the other phases of ice due to its close structural similarity to sI hydrate at the molecular level: the average hydrogen bond length in the hydrate is only 1% longer on average than in ice I h and the hydrate O-O-O angles differ from the tetrahedral angles of ice I h by only 3.7 • . 1 In the same manner that we calculated ∆E sI coh (a 0 ) for sI methane hydrate, we have also calculated ∆E ice coh for our test set of xc functionals and FFs by fitting the cohesive energy of the bulk ice I h crystal to Murnaghan's equation of state (a more comprehensive overview of the ice results is given in the SI). For DMC, we use the value of ∆E ice coh reported in Ref. 16. From these calculations, we also obtain the equilibrium volume of the ice I h crystal. In Fig. 4 we show the difference in computed volume using DFT/FF from that using DMC for sI methane hydrate, plotted against the same quantity for ice I h . There is a strong correlation between the errors in computed volumes for the hydrate and ice I h , suggesting that the primary factor in obtaining reasonable lattice volumes for the hydrate is an accurate description of the hydrogen bonded water framework. We also show in Fig. 4 the differences in DFT/FF values for ∆E sI coh and ∆E empty coh from DMC, again plotted against the DFT/FF-DMC difference for the ice I h cohesive energy. As for the volumes, we see a very strong positive correlation between the errors in the hydrate cohesive energies and those for ice I h . In fact, there is a near perfect correlation for the error ∆E empty coh and the error in ∆E ice coh , the significance of which will become apparent when we look at the dissociation behaviour of the hydrate to ice I h and methane gas.
In comparing the cohesive energies of the DFT xc functionals and point charge FFs to DMC, we are taking the vapour phase of both methane and water as our reference state. More important to the phase equilibria of gas hydrates, however, is the relative energy of the hydrate with respect to methane gas and another condensed phase of water, either liquid or ice. 1 Whilst the cancella- tion of errors in ∆E CH 4 and ∆E empty coh means that PBE has a good overall agreement with the DMC cohesive energy, it is straightforward to demonstrate that the error in ∆E CH 4 arising from the neglect of vdW interactions can lead to severe consequences regarding the thermodynamic stability of sI methane hydrate. Consider the process of sI methane hydrate dissociating to ice I h and methane gas: The associated energy cost ∆E sI→ice diss can be computed as (see SI): The results of these calculations are presented in Table I. It is clear that the results for PBE are disastrous: sI methane hydrate is unstable with respect to dissociation to ice I h and methane gas by 67 meV/CH 4 (i.e. it is 67 meV/CH 4 exothermic). In contrast, DMC predicts dissociation to be an endothermic process, costing 155 ± 34 meV/CH 4 . We note here that the experimental enthalpy of dissociation 68 at standard temperature and pressure is 188 ± 3 meV/CH 4 suggesting that the DMC value is reasonable. 69 All of the dispersion-corrected functionals improve on the GGA functional in this respect, predicting that the hydrate is stable with respect to ice and methane gas. PBE-D2 gives the best agreement with DMC, followed by LDA and PBE-vdW TS , although it should be kept in mind that these calculations have been performed at the equilibrium volume of the xc functional used. The vdW-DFs over-stabilise the hydrate by a factor of approximately two. Unsurprisingly, the trends in ∆E ice coh closely follow those of ∆E CH 4 , with the errors in describing the hydrogen bonded water network more-orless cancelling between ∆E empty coh and ∆E ice coh (as shown in Fig. 4). As such, the FFs also give good agreement with the DMC result. The fact that the point charge FFs predict ∆E sI coh to be too exothermic can be attributed to the enhanced dipole moment of the isolated water molecules in these types of potentials, 70,71 which has been shown to lead to too high vaporisation enthalpies of ice I h for the TIP4P-2005 and TIP4P-ICE potentials. 72 Indeed, Vega and co-workers 72 have found that it is impossible to simulataneously fit the melting temperature of ice I h and the enthalpy of vaporisation for such models. It is therefore probably expecting too much of the rigid point charge FFs to give reasonable results for both ∆E sI coh and ∆E sI→ice diss whilst also maintaining favourable densities and coexistence/melting temperatures for the hydrate and ice I h . 73 Use of an explicitly polarizable water potential may go some way to improving this situation. 74

IV. CONCLUSIONS
We have presented high-quality DMC reference data for bulk sI methane hydrate and evaluated the performance of several commonly used xc functionals and point charge force fields. We have found that none of the DFT methods tested give particularly satisfactory results. We find that vdW forces are crucial to the stability of methane hydrate with respect to dissociation to ice I h and methane gas, although the vdW-DF flavour of xc functionals over-stabilise the hydrate by approximately a factor of two. This effect is less severe with the PBE-D2 and PBE-vdW TS functionals, although their equilibrium volumes are too small compared to DMC. PBE, which neglects dispersion interactions, incorrectly predicts that methane hydrate is unstable with respect to dissociation to ice and methane gas. By overbinding the hydrogen bonded water framework, PBE's poor description of the water-methane interaction is compensated, giving a good overall agreement with the DMC cohesive energy of the bulk hydrate. This last point highlights the difficulty that DFT xc functionals face in describing mixed phase systems such as gas hydrates; in order to obtain a good overall description, it is necessary to be able to accurately describe both the hydrogen bonded water framework and the dispersion bound methane. We have also seen that point-charge, all-atom force fields tend to overbind the hydrate lattice, although their agreement with DMC for the dissociation energy to ice and vapour, and for the structure for the bulk crystal, is good. From our knowledge of the literature 72 on the performance of simple point charge FFs for ice, it is unlikely that such FFs will ever simultaneously describe both the cohesive energy of the hydrate crystal and the energetics of dissociation to other condensed phase water systems.
Earlier in this article, we remarked that the 23:4 ratio of water to methane amplified the apparently small errors in the water-water interactions compared to the water-methane interactions. We also saw that the high water content means that the errors in describing the hydrate are strongly correlated to the errors in describing ice I h . However, such a high water-methane ratio also means that there is the possibility for significant many body interactions between the methane and water (e.g. a single isolated 5 12 cage has 190 water-methane-water triplets). Indeed, a separate independent study investigating the binding energy of methane to a gas phase 5 12 cage through a many-body expansion of the total energy has found significant contributions to the DFT error beyond those in the two-body interactions (symmetry adapted perturbation theory calculations also showed that the DFT methods have insufficiencies other than those associated with the neglect of long-range dispersion interactions). 75 Although we have not considered the effects of including exact exchange, it is unlikely that this will significantly improve the DFT description of sI methane hydrate. For example, using the PBE0 76 results for ice I h from Refs. 17, and for a methane molecule binding to a gas phase 5 12 water cage from Ref. 75, we predict that this hybrid xc functional will give a reasonable prediction of the hydrate structure (similar to PBE) but will still incorrectly destabilise the hydrate with respect to methane gas and ice I h . Including dispersion corrections to this functional, such as PBE0-D2 or PBE0-vdW TS , can therefore be expected to also give similar results to PBE-D2 and PBE-vdW TS . It has recently been shown 77 that accurate DMC reference data in combination with Gaussian approximation potentials 78 can be used to systematically correct the 'beyond two-body' errors associated with GGA functionals for water nano-droplets and bulk liquid water. Such an approach is also likely to be more successful in improving the performance of DFT xc functionals for gas hydrates compared to the pairwise additive dispersion corrections examined here.