Jérôme Cuny

Franck Jolibois

Iann C Gerber
email: igerber@insa-toulouse.fr

Evaluation of Gas-to-Liquid 17 O Chemical Shift of Water: a Test Case for Molecular and Periodic Approaches

Modelling liquid water features is a challenging and ongoing task that brings together a number of computational issues related to both the description of its electronic and geometrical structure. In order to go a step further in the understanding of this peculiar liquid, we present a thorough analysis of NMR gas-to-liquid 17 O and 1 H shifts of water using density functional theory based molecular dynamics. In order to be as consistent as possible, we consider the influence of basis sets, exchangecorrelation functionals and structural models, in both molecular and periodic schemes to evaluate 17 O and 1 H nuclear shieldings. We show that strong error compensations between functional and basis-set expansion can be obtained in molecular approaches which artificially produces good 17 O gas-to-liquid shifts with relatively small basis sets.

We also demonstrate that despite their ability to provide reliable liquid phase structures, generalized-gradient approximation based exchange-correlation functionals lead 1 to strongly inconsistent values for 17 O gas-to-liquid shift. This latter property is shown to be strongly influenced by intra-molecular electronic delocalization, accentuated by the surrounded molecules. In contrast, 1 H is less sensitive to this effect. By including a Hartree-Fock exchange term, through the use of hybrid functionals which partially correct the self-interaction error, better agreement with experimental values is obtained.

The present study provides a detailed guideline to properly evaluate gas-to-liquid shifts in hydrogen bonded systems and emphasizes that, for nuclear shieldings, an accurate electronic structure evaluation prevails over the description of the liquid structure.

Introduction

Nuclear magnetic resonance (NMR) is a key spectroscopic method to characterize the structure and dynamics of molecular liquids, solutions or solids. [START_REF] Emsley | Milestones in the First fifty Years of NMR[END_REF][START_REF] Emsley | Forty years of Progress in Nuclear Magnetic Resonance Spectroscopy[END_REF] However, despite the great achievements in the development and optimization of this tool, an unambiguous analysis of the experimental spectra is often difficult to achieve. In spite of the great sensitivity of NMR to local environment, this can lead to a great loss of chemical information. To overcome this limitation, NMR has been efficiently complemented by quantum chemical calculations. Such approaches have emerged a long time ago with the pioneering studies by Lipscomb and coworkers, [START_REF] Stevens | Perturbed Hartree-Fock Calculations. I. Magnetic Susceptibility and Shielding in the LiH Molecule[END_REF] and have been developed in both molecular finite size and periodic frameworks. In the former case, the main theoretical limitation lies in the definition of the gauge origin. A naive and ill-suited definition can lead to major convergence issues with respect to the basis set. This inevitably limits both the accuracy of the results and the system sizes that can be handled. Nowadays, it is well recognized that the Gauge-Including Atomic Orbitals (GIAO) formalism partly solves this question. [START_REF] Ditchfield | Molecular Orbital Theory of Magnetic Shielding and Magnetic Susceptibility[END_REF][START_REF] Wolinski | Efficient Implementation of the Gauge-Independent Atomic Orbital Method for NMR Chemical Shift Calculations[END_REF] For systems described within periodic boundary conditions, other difficulties arise from the combination of both the infinite nature of the systems and the plane-wave pseudopotential scheme usually applied to compute electronic structures. In the early 2000's, Pickard and Mauri proposed the Gauge-Including Projector Augmented Wave (GIPAW) approach to tackle this issue. [START_REF] Pickard | All-Electron Magnetic Response with Pseudopotentials: NMR Chemical Shifts[END_REF] For more details on the GIPAW approach and its applications, see Ref.

7, 8 and references therein.

Both computational means provide isotropic nuclear shieldings (σ iso ) that cannot be directly compared to experimental chemical shifts since the latter are defined as a difference between nuclear shieldings of a probed nucleus and the one of the same nucleus in a reference compound. As a consequence, computational schemes can only provide relative chemical shifts as long as the nuclear shielding of the reference compound is not properly evaluated. In [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O NMR experiments, the reference system is usually liquid water and an accurate theoretical evaluation of its nuclear shielding is thus desirable. The ability of a given methodology in determining this quantity can be probed by comparing theoretical gas-to-liquid shifts of both [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H isotropic nuclear shielding of water, defined as ∆σ iso theo = σ iso liquid -σ iso gas ,

to experimental values. For 1 H and [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O, ∆σ iso exp are -4.3 and -36.1 ppm, respectively. [START_REF] Hindman | Proton Resonance Shift of Water in the Gas and Liquid States[END_REF][START_REF] Florin | 17 O NMR Shifts in H 17 2 O Liquid and Vapor[END_REF][START_REF] Raynes | The 17 O Nuclear Magnetic Shielding in H 2 17 O and D 2 17 O[END_REF][START_REF] Modig | Proton Magnetic Shielding Tensor in Liquid Water[END_REF][START_REF] Wasylishen | A Revised Experimental Absolute Magnetic Shielding Scale for Oxygen[END_REF] However, this is not an easy task since nuclear shielding evaluation issues are supplemented by the difficulty to properly describe the structure of liquid water. At the density functional theory (DFT) level, this latter question has been an ongoing debate for more than two decades. Indeed, since the pioneering work by Laasonen et al. in the mid-1990's, [START_REF] Laasonen | Ab initio" Liquid Water[END_REF] a number of studies have been performed to rationalize the strengths and weaknesses of DFT for the structural and dynamical description of liquid water in ab initio molecular dynamics (MD). Obtaining a reasonable structure is a necessary requirement to further model the spectroscopic features of liquid water at the DFT level of theory, in particular NMR properties.

In this context, the ability of a given computational methodology to reproduce the 17 O (and 1 H) ∆σ iso exp of water represents a natural test-case for its accuracy as ∆σ iso is directly related to the nature and strength of intermolecular interactions within the liquid. Hence, it is a real fingerprint on how good the description of the chemical and electronic environment of a water molecule in the liquid is. Unfortunately, this means that ∆σ iso is extremely sensitive to various computational parameters: basis-set definition and size, exchange-correlation (xc)-functional, accuracy of the structural model and its statistical convergence.

Various theoretical studies have already been conducted to evaluate [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H ∆σ iso in water. Using molecular formalisms, water clusters of various sizes were extracted from MD simulations to explicitly account for solvent effects, [START_REF] Malkin | Solvent Effect on the NMR Chemical Shieldings in Water Calculated by a Combination of Molecular Dynamics and Density Functional Theory[END_REF][START_REF] Cybulski | On the Calculations of the Nuclear Shielding Constants in Small Water Clusters[END_REF] when implicit models were tested too, [START_REF] Klein | Ab Initio Calculations of O-17 NMR-Chemical Shifts for Water. The Limits of PCM Theory and the Role of Hydrogen-Bond Geometry and Cooperativity[END_REF] to evaluate nuclear shieldings at various levels of theory going from Hartree-Fock to DFT. [START_REF] Pennanen | Nuclear Magnetic Resonance Chemical Shifts and Quadrupole Couplings for Different Hydrogen-Bonding Cases Occurring in Liquid Water: A Computational Study[END_REF] The influence of the xc-functional was reported. [START_REF] Fileti | Isotropic and Anisotropic NMR Chemical Shifts in Liquid Water: A Sequential QM/MM Study[END_REF] Overall, accurate results can only be yielded by taking advantages of error's cancellation originating from (i) the type of potential used for the MD simulation, (ii) the size of the selected clusters used for the NMR parameter calculations, (iii) the basis-set extension, and finally (iv) the choice of the xc-functional.

Further studies were also performed using periodic boundary conditions. The first one was conducted by Pfrommer et al. [START_REF] Pfrommer | NMR Chemical Shifts of Ice and Liquid Water: The Effects of Condensation[END_REF] using the so-called Mauri, Pfrommer and Louie (MPL)

formalism. [START_REF] Mauri | Ab Initio Theory of NMR Chemical Shifts in Solids and Liquids[END_REF] In this study, a 300 K DFT-based MD simulation using the BLYP [START_REF] Becke | Density-Functional Exchange-Energy Approximation With Correct Asymptotic Behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density[END_REF] xcfunctional was used to generate the structural model. Nuclear shieldings were subsequently evaluated using the local density approximation. In this study, the error induced by the pseudo-potential approximation on the valence wave-functions in the core regions was considered constant from gas to liquid. In spite of this, the author calculated [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H ∆σ iso theo in quantitative agreement with experimental values, without considering any statistical approach. Other studies were further conducted although they mainly focused on the 1 H nuclear shielding of water under various conditions. [START_REF] Ceriotti | Nuclear Quantum Effects and Hydrogen Bond Fluctuations in Water[END_REF][START_REF] Sebastiani | Ab-initio Study of NMR Chemical Shifts of Water Under Normal and Supercritical Conditions[END_REF][START_REF] Murakhtina | First Principles and Experimental 1H NMR Signatures of Solvated Ions: The Case of HCl(aq)[END_REF][START_REF] Schmidt | Bulk Chemical Shifts in Hydrogen-Bonded Systems from First-Principles Calculations and Solid-State-NMR[END_REF][START_REF] Banyai | NMR Chemical Shifts as a Tool to Analyze First Principles Molecular Dynamics Simulations in Condensed Phases: the Case of Liquid Water[END_REF] Overall, although those studies have provided a number of important achievements, to the best of our knowledge, none of them provided a consistent picture about the interplay between the computational parameters influencing the [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H gas-to-liquid shifts of water.

In the present study, we intend to provide a comprehensive view on how to evaluate the [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H NMR gas-to-liquid shifts of water using DFT formalism only. As error compensations can lead to acceptable but non-physical results, a proper evaluation of the impact of various computational parameters: basis-set extension, xc-functional approximation, amount of Hartree-Fock (HF) exchange in the case of hybrid xc-functional, [START_REF] Becke | A New Mixing of Hartree-Fock and Local Density-Functional Theories[END_REF] GIAO in conjunction with a cluster ansatz vs fully periodic GIPAW approach, and their mutual interplay is provided. This allows us to draw a reliable picture on how correct [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H gas-to-liquid shifts of water can be obtained and be served as reference values for further studies. The outline of the article is as follows: the computational methods employed along the article are described in Section and the results are presented and discussed in Section .

This section is divided into three parts: description and validation of the molecular dynamics simulations, and then NMR results obtained from molecular and periodic calculations are detailed. The main outcomes and perspectives are summarized in the conclusion.

Computational Details

Periodic MD Simulations. We carried out periodic DFT simulations of liquid water using the Vienna ab initio simulation package (VASP). [START_REF] Kresse | Ab Initio Molecular Dynamics for Liquid Metals[END_REF][START_REF] Kresse | Efficient Iterative Schemes for Ab Initio Total-Energy Calculations using a Plane-Wave basis Set[END_REF] The code uses the full-potential

Projector Augmented Waves (PAW) framework. [START_REF] Blöchl | Projector Augmented-Wave Method[END_REF][START_REF] Kresse | From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method[END_REF] Standard versions of the PAW potential for both H and O atoms were used in conjunction with an energy cut-off of 400 eV. Sampling of the Brillouin zone was reduced to Γ-point only. Exchange-correlation interactions were approximated by using either the non-spin polarized version of the Generalized Gradient Approximation (GGA) PBE, [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] or the optB88-vdW scheme [START_REF] Dion | Van der Waals Density Functional for General Geometries[END_REF][START_REF] Román-Pérez | Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes[END_REF][START_REF] Klimeš | Chemical Accuracy for the van der Waals Density Functional[END_REF][START_REF] Klimeš | Van der Waals Density Functionals Applied to Solids[END_REF] in order to evaluate the impact of weak intermolecular interactions on the liquid structure, as suggested in Ref.

29,39. Born-Oppenheimer MD simulations were carried out under constant volume and temperature using a Nosé thermostat defined by a period of 20 fs. [START_REF] Nosé | A Unified Formulation of the Constant Temperature Molecular Dynamics Methods[END_REF] The simulations consisted in 64 water molecules in a cubic unit-cell of 12.42 Å which corresponds to a density of 1 g.cm -3 . Equations of motion were integrated with a 0.5 fs time step for a total of 20 ps for each simulation. Data were extracted after an equilibration time of 10 ps. By using this time scale, a reasonable agreement with longer simulations is obtained in terms of positions and heights of the extrema of the pair radial distribution functions. [START_REF] Pham | Structure and Dynamics of Aqueous Solutions from PBE-Based First-Principles Molecular Dynamics Simulations[END_REF] Simulations were performed at 300 and 330 K, the latter being used to artificially mimic the fluctuations resulting from nuclear quantum effects (NQEs) as previously proposed. [START_REF] Morrone | Nuclear Quantum Effects in Water[END_REF] In the following, three MD trajectories denoted MD-PBE-300, MD-OptB88-300 and MD-PBE-330 were used as basis for NMR calculations. To do so, 50 randomly selected frames were considered from each MD simulation for subsequent GIPAW and GIAO NMR calculations.

Periodic NMR Calculations. Periodic calculations of 1 H and 17 O nuclear shieldings were performed within the GIPAW approach, [START_REF] Pickard | All-Electron Magnetic Response with Pseudopotentials: NMR Chemical Shifts[END_REF][START_REF] Yates | Calculation of NMR Chemical Shifts for Extended Systems using Ultrasoft Pseudopotentials[END_REF] as implemented in the VASP package. [START_REF] Vasconcelos | Finite-Field Implementation of NMR Chemical Shieldings for Molecules: Direct and Converse Gauge-Including Projector-Augmented-Wave Methods[END_REF] Several xc-functionals were used for those calculations: PBE, optB88-vdW, B3LYP, 72 PBE0, [START_REF] Perdew | Rationale for Mixing Exact Exchange with Density Functional Approximations[END_REF][START_REF] Adamo | Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0 Model[END_REF] and a modified version of PBE0 which includes 50% of Hartree-Fock (HF) exchange to investigate the influence of the amount of exact exchange in the xc-term. Indeed, it has been recently shown in the case of [START_REF] Ramirez | Quantum Corrections to Classical Time-Correlation Functions: Hydrogen Bonding and Anharmonic Floppy Modes[END_REF] F, that this parameter is a critical issue in the calculation of nuclear shieldings. [START_REF] Alkan | Role of Exact Exchange and Relativistic Approximations in Calculating 19 F Magnetic Shielding in Solids Using a Cluster Ansatz[END_REF] The HSE xc-functional was also tested, [START_REF] Heyd | Assessment and Validation of a Screened Coulomb Hybrid Density Functional[END_REF][START_REF] Heyd | Energy Band Gaps and Lattice Parameters Evaluated with the Heyd-Scuseria-Ernzerhof Screened Hybrid Functional[END_REF][START_REF] Paier | Screened Hybrid Density Functionals Applied to Solids[END_REF] since it is a range separated hybrid-type xc-functional that includes a 25% short range HF exchange term that allows to discriminate between short and long range contributions to nuclear shieldings.

Convergence issues with respect to the extension of the plane-wave basis set were carefully tested. The results are displayed in the left panels of Figure 1 for an isolated water molecule using PBE and PBE0 xc-functionals. The right panels of Figure 1 All GIPAW calculations were performed using a 1 × 1 × 1 k -point grid for the integration of the Brillouin zone. The use of denser grids has a very limited impact on the result's accuracy. Indeed, considering a single 64-molecules snapshot, average variations of less than 0.1 and 2.5 ppm for 1 H and 17 O nuclear shieldings, respectively, are observed when increasing the size of the grid from 1 × 1 × 1 to 4 × 4 × 4. Harder PAW potentials were also tested on the same snapshot with an energy cut-off of 800 eV leading to upshifts of only 0.9 and 0.7 ppm for 1 H and 17 O nuclear shieldings, respectively. Molecular NMR Calculations. Nuclear shieldings were computed on molecular clusters using the GIAO approach 4,5 implemented in the Gaussian software package. 79 Four among the six xc-functionals (PBE, PBE0, HSE and B3LYP) used in periodic calculations were employed together with several Popple-type basis sets with increasing size to investigate basis-set dependency. It is worth pointing out that it has been largely reported that NMR parameter accuracy are strongly basis-set dependent. [START_REF] Jensen | Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods[END_REF][START_REF] Kupka | Convergence of Nuclear Magnetic Shieldings in the Kohn-Sham Limit for Several Small Molecules[END_REF][START_REF] Flaig | Benchmarking Hydrogen and Carbon NMR Chemical Shifts at HF, DFT, and MP2 Levels[END_REF] The convergence of the calculated nuclear shieldings with respect to the size of the con-sidered cluster was first investigated. We extracted, from the MD-PBE-300 run, 1025 H 2 O-(H 2 O) n distinct clusters with n going from 1 to 36. Those 1025 structures were obtained by considering 25 different central water molecules along 41 distinct frames. The corresponding σ iso O and σ iso H obtained at the B3LYP/6-311G(d,p) level are displayed in Figure 2. σ iso H is fully converged for n = 24, when σ iso O appears more difficult to converge as a small but continuous decrease in the average value is observed even for n greater than 24. However, between n = 16 to 36, the variation is only 2.4 ppm which is small as compared to the For the same set of computational parameters, i.e. B3LYP/6-311G(d,p), we also tested the nuclear shielding convergence with respect to the number of considered configurations.

To do so, we calculated the mean, standard deviation and standard deviation over the mean obtained when increasing the total number of frames considered in the statistical averaging. The standard deviation over the mean as to be understood as follows: at point k on the x-axis, k mean values are calculated (for 1 selected configuration only to k selected configuration), a standard deviation can then be calculated from those k mean values. This provides an information about how the mean value fluctuates when an increasing number of configurations is considered. In contrast to the standard deviation which converges towards a non-zero value, the standard deviation over the mean tends to 0 when an infinite number of configurations is considered in the averaging. It is thus a better indicator of the convergence of the statistical sampling. The results are reported in Figure 3 for the same 1025 distinct

H 2 O-(H 2 O) 36 clusters discussed above.
As a consequence, the mean and standard deviation values obtained for 1025 selected configurations are equal to the values reported in Figure 2 for n = 36. Although the standard deviation is quite large, highlighting the variety of chemical environments existing around a given kind of atom, the standard deviation over the mean converges rapidly for both atoms. When all the frames are considered, it is equal to 0.17 xc-Functional and Basis-Set Issues in GIAO Calculations. As the choice of xcfunctional and basis-set is crucial in NMR GIAO calculations we performed various tests to check their respective influence on the calculated 1 H and 17 O nuclear shieldings. In that case, to avoid performing statistical sampling on a large series of water clusters with a number of different computational settings, basis-set and xc-functional influence on nuclear shieldings were evaluated on a single H 2 O-(H 2 O) 42 water cluster randomly extracted from the MD-PBE-300 simulation. A large enough cluster was considered in order to avoid any finite-size effect. In the Supporting Information, calculations performed on eleven configurations and leading to average 1 H and [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shieldings are also presented (see Figure S2 andS3). They lead to very similar results as the ones discussed below obtained from the single 311G(d,p) to 6-311++G(2d,2p), PBE, B3LYP, HSE and PBE0 always display an increasing agreement with experiment. In contrast, for 6-31G(d,p) and 6-31G(2d,2p), PBE leads to the best results while B3LYP, HSE and PBE0 perform equally well. Considering the largest basis set, 6-311++G(2d,2p), the PBE result is equal to -62.2 ppm. This value is very close to the statistically converged GIPAW results obtained with the same xc-functional as discussed in section . This suggests that both approaches converge to similar 17 O gas-to-liquid shift despite their differences. Finally, to ensure that the present results do not depend on the basis set construction, additional basis sets were tested. The results are presented in the basis sets lead to similar results.

H 2 O-(H 2 O) 42 configuration.
For 1 H, differences lower than 2 ppm as compared to experimental data can only be obtained with the 6-31G(d,p) basis set, see Figure 4. All other basis sets show errors greater than 2 ppm which corresponds to a ∼50 % of difference as compared to experiment. In contrast to [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O, Figure 4 also demonstrates that the xc-functional influence is rather limited for 1 H, although the inclusion of a HF exchange term systematically slightly improves the results. Indeed, PBE always leads to the worst results, while B3LYP, HSE and PBE0 perform equally well. In the same way, diffuse functions have a very limited impact on the 1 H gasto-liquid shift. Despite displaying a smaller range, 1 H nuclear shielding also appears to be converged for the largest basis sets: 6-311+G(d,p), 6-311++G(d,p) and 6-311++G(2d,2p), and the corresponding value at the PBE level is close to the GIPAW results (see discussion in section ).

In the light of those preliminary results performed on a single water cluster, we decided to compute all subsequent GIAO nuclear shieldings using B3LYP. B3LYP was chosen as it provides 1 H and 17 O gas-to-liquid shifts comparable to HSE and PBE0 and is widely used in the quantum chemistry community. In terms of basis sets, we considered two small basis sets, 6-31G(d,p) and 6-31G(2d,2p), and two larger basis sets: 6-311G(d,p) and 6-311++G(d,p), in order to probe their respective behavior with respect to statistical sampling.

Results and Discussion

MD Simulations

As stated in introduction, a first and important step of the present study is to provide a correct structure of liquid water for subsequent evaluation of NMR parameters. However, this is not an easy task as an accurate ab-initio modeling of liquid water is difficult to conduct. Indeed, despite the large number of studies dedicated to this question, 15-44 few of them have been able to succeed in a quantitative way. One of the best achievements was obtained by DiStasio et al. [START_REF] Distasio | The Individual and Collective Effects of Exact Exchange and Dispersion Interactions on the Ab Initio Structure of Liquid Water[END_REF] who conducted a simulation of 64 water molecules at 330 K at the PBE0+TS-vdW(SC) level, i.e. using a hybrid PBE0 xc-functional, [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF][START_REF] Adamo | Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Models[END_REF] in conjunction with a self-consistent dispersion-correction term. [START_REF] Tkatchenko | Accurate Molecular Van Der Waals Interactions from (85)[END_REF] The authors obtained oxygenoxygen (g OO (r)) and oxygen-hydrogen (g OH (r)) radial distribution functions as well as O-O-O angular distribution functions that are in very good agreement with the experimental data of Soper and Benmore. 85 In the present work, we do not intend to reproduce those calculations that are computationally extremely expensive. Instead, we resort to more traditional xcfunctionals: PBE and optB88-vdW, the latter being used to produce a MD trajectory which suffers less of the famous liquid water over-structuration issue resulting from the use of standard GGA xc-functionals. [START_REF] Sit | Static and Dynamical Properties of Heavy Water at Ambient Conditions from First-Principles Molecular Dynamics[END_REF] We also present the results of a PBE-run performed at 330 K (MD-PBE-330).

The g OO (r) and g OH (r) radial distribution functions for the three simulations are presented in Figure 5. As expected, the PBE xc-functional at 300 K leads to overstructured radial distribution functions, for both oxygen-oxygen and oxygen-hydrogen pairs, in agreement with the literature. [START_REF] Sit | Static and Dynamical Properties of Heavy Water at Ambient Conditions from First-Principles Molecular Dynamics[END_REF] The MD-PBE-330 and MD-OptB88-300 runs produce very similar curves which are significantly softened when compared to the MD-PBE-300 ones, i.e. closer to the experimental results. [START_REF] Soper | The Radial Distribution Functions of Water as Derived from Radiation Total Scattering Experiments: Is There Anything We Can Say for Sure[END_REF] Indeed, for g OO (r), the intensity of the first peak is reduced by approximately 0.6 (MD-PBE-330) and 0.7 (MD-OptB88-300), while the intensity of the first minimum is increased by 0.17 (MD-PBE-330) and 0.24 (MD-OptB88-300), respectively. An increase in the positions of the first maximum and first minimum, by 0.04 Å with respect to the MD-PBE-300 trajectory is also observed for both MD-OptB88-300 and MD-PBE-330.

One can thus expect that such improvements of the general structure of the simulated liquid water may influence NMR parameters.

Cluster Approximation Calculations

Reference Water Molecule. To discuss the gas-to-liquid shift of water, we need to evaluate both the nuclear shielding of the water molecule in the gas and liquid phases. In the former case, two approaches can be used: (i) consider a water molecule at equilibrium geometry for the considered xc-functional/basis-set combination; (ii) perform a statistical average over a MD run on a single water molecule in vacuum. The second approach advantage is that it takes into account the vibrational effects in the gas phase. However, those effects are expected to be rather small due to the strong harmonic character of the vibrational modes of the isolated water molecule. As shown in the Supporting Information, the difference between the two approaches is negligible: 0.08 and 1.8 ppm for σ iso H and σ iso O , respectively. Consequently, in the following, we consider the nuclear shielding of an isolated water molecule in its equilibrium geometry as the gas-phase reference for both GIPAW and GIAO calculations.

Molecular Dynamics Simulation and Basis-Set Influence. [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] basis sets provide ∆σ iso theo values that are too small in absolute value for [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O, and too large for 1 H, while 6-311G(d,p) and 6-311++G(d,p) leads to values that are too large for both [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H. This is fully in line with the results of Figure 4 obtained on a unique water cluster. Indeed, for 1 H, both small and large basis sets lead to similar results, the maximum difference is observed between 6-31G(2d,2p) and 6-311G(d,p). For 17 O, 6-31G(2d,2p)

provide results that are closer to the experimental value while 6-311++G(d,p), the largest basis set considered here, lead to the largest difference. Consequently, statistical averaging confirms the drastic influence of the basis-set on the [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shielding regardless of the structural model and its more limited impact on the 1 H nuclear shielding.

Considering the MD settings and the 6-31G(d,p) and 6-31G(2d,2p) basis sets, when the temperature is raised from 300 to 330 K using the PBE xc-functional, the theoretical error decreases by more than 4.0 and 1.5 ppm for oxygen and hydrogen, respectively. For 17 O, no significant difference appears between the two simulations (MD-PBE-300 and MD-OptB88-300) performed at 300 K. This is surprising as MD-OptB88-300, which includes van der Waals corrections, leads to g(r) functions closer to experiments. The same observation holds for 1 H, although larger differences than for this seeming agreement between theory and experiment is to be taken with care seeing the limited convergence of those two basis sets (see section .)

In conclusion, as observed in the previous section, GIAO results for [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shieldings are extremely sensitive to basis set effects. Indeed, large basis sets lead to the largest discrepancy with experimental data whereas small basis sets lead to values surprisingly close to the experimental value. This fortuitous result will be further analysed in section in the light of GIPAW calculations. In addition, the present results suggest that difference between the structural models are limited if compared to the impact of the basis set although an increase of the temperature by only 30 K appears to noticeably modify the 1 H nuclear shieldings. In these circumstances, extrapolation and/or rationalization of GIAO results are almost impossible and appear to be highly system-dependent. Periodic Calculations.

Statistical Convergence. One advantage of periodic over molecular calculations is that finite-size effects are intrinsically avoided and large statistical approach can be assessed, since GIPAW calculations over 50 frames lead to 6400 and 3200 1 H and [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shielding values, respectively. This is three times larger than for the molecular cluster approach discussed above. One can expect that the present nuclear shielding calculations are fully converged in terms of considered structures' number. To further support this assertion, Structural Model and Density Functional. Table 1 contains the mean gas-to-liquid shift values obtained from twelve different set of computational parameters: two different MD simulations (MD-OptB88-300 and MD-PBE-330), and for each one, six different xc-functionals used for GIPAW calculations on 50 different frames. In Table S1 of the Supporting Information, data obtained from the MD-PBE-300 trajectory are also provided.

Interestingly the effect of the structural model is very limited. Indeed, if one compares MD-PBE-300 with MD-PBE-330 results, increasing the temperature only leads to a small decrease of all nuclear shielding values, with largest absolute variations of 2.3 and 0.5 ppm for [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H, respectively, in the case of the PBE xc-functional. In the same way, van der Waals corrections applied to the MD simulations do not significantly affect final nuclear shielding values, when it impacts positively the liquid structure. Indeed, the largest differences between MD-OptB88-300 and MD-PBE-300 are obtained with PBE and are equal to 0.7 and 3.4 ppm only (see Table 1 and Table S1) for 1 H and 17 O, respectively. Besides, the largest differences between MD-OptB88-300 and MD-PBE-330 are only 0.2 and 1.1 ppm for 1 H and 17 O. Those results demonstrate that MD simulations leading to different pair radial distribution functions (for instance MD-OptB88-300 and MD-PBE-300) lead to very similar gas-to-liquid shifts for both 1 H and 17 O when considering the same xc-functional for GIPAW calculation. The same conclusion also holds for MD simulations leading to very similar pair radial distribution functions as for instance MD-OptB88-300 and MD-PBE-330. This statement is consistent with the results of section when considering the GIAO molecular approach.

The downside is that the choice of xc-functional for GIPAW calculations is crucial. For the MD-OptB88-300 trajectory, the variation amplitude of the 17 O gas-to-liquid shift is large when considering GGA xc-functionals and hybrid ones: from -64.7 ppm for PBE to -23.1 ppm for PBE0-50%. A more detailed analysis shows that such a large variation originates from both terms involved in ∆σ iso theo . Indeed, in Table S2 and S3 of the Supporting Information, we report the σ iso gas values used to calculate the gas-to-liquid shifts provided in Table 1. The amplitude of σ iso gas values for [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shielding is of only ∼20 ppm, from 320.9 ppm for PBE to 340.6 ppm for HSE xc-functional. As a consequence, the interaction between the water molecules has a strong impact on the calculated [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shieldings and the six considered xc-functionals do not provide an equivalent description of this interaction. This demonstrates that the electronic density in the close vicinity of the nucleus of interest has a strong impact on its nuclear shielding. It also shows that by including an HF exchange term, the hydrogen shielding decreases while the oxygen shielding increases.

This effect is directly related to the well-known tendency of GGA xc-functionals to overdelocalize the electronic density, resulting from intrinsic self-interaction error, 88 when hybrid xc-functionals partially correct it. Interestingly, GGA xc-functionals provide smaller [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shieldings as compared to hybrid ones even for a single water molecule which further highlights this tendency. The situation is inverted for 1 H. To illustrate this point, in the top panel of Figure 8, the difference between the PBE0 and PBE electronic densities of an isolated water molecule is depicted. The re-localization of the density on the oxygen lone pairs is clearly demonstrated. This is in line with the very recent finding of Ref. 44 which shows that the center of oxygen lone pairs in water tend to be closer to the nuclei when a meta-GGA xc-functional is used to describe exchange-correlation terms. This intramolecular electronic density re-organization resulting from the inclusion of HF exchange term is even more pronounced when the liquid phase is considered, see bottom panel of Figure 8, since hydrogen bonds are formed with a strong interaction between the oxygen lone pairs and the first neighbouring hydrogen atom.

We show here that an accurate description of the hydrogen bonds is of paramount importance for nuclear shielding calculation in such hydrogen-bonded systems. From the experiments. This indicates that the balance between HF and DFT-based exchange term is subtle. When too much HF exchange is taken into account, the re-localization of the density is too pronounced, making the oxygen σ iso liquid value too close to the gas phase.

This artificially limit the impact of the the chemical environment on the nuclear shielding of oxygen. Mind that this over-relocalization has its counterpart on the hydrogen nuclear shielding, the shielding being weaker for PBE0-50%. HSE and B3LYP seems to be more balanced and lead to similar values slightly below the experimental one. The best agreement for [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O is obtained with the PBE0 xc-functional, in combination with the MD-OptB88-300 simulation, that leads to a ∆σ theo -∆σ exp value equals to -1.5 ppm only. This is slightly better than the results obtained from the GIAO-cluster-based approach discussed in section and obtained with small basis sets. However, GIPAW is not subject to basis-set effects as for the GIAO calculations and the comparison of the two approaches allow to partly rationalize them. Indeed, better agreement with experiment is obtained with the smallest basis sets at the PBE level as they do not provide enough flexibility to the electronic density to delocalize.

So, this artificially leads to results close to experiment. When large basis sets are considered, in particular when diffuse functions are added on the oxygen atoms, description of electronic Considering the 1 H gas-to-liquid shift, the situation is slightly different. The ∆σ iso H varies from -7.0 (PBE) to -5.4 ppm (PBE0-50%), meaning that when more HF exchange is included, better the agreement with the experimental value is. However, this latter is not reached even with the PBE0-50% xc-functional which suggests that some ingredients could be missing in the present simulations. In particular, Ceriotti et al. demonstrated that NQEs can change computed 1 H chemical shifts by as much as 0.5 ppm. [START_REF] Ceriotti | Nuclear Quantum Effects and Hydrogen Bond Fluctuations in Water[END_REF] Consequently, a subtle balance between the inclusion of HF exchange and the description of NQEs may be necessary to exactly reproduce the 1 H gas-to-liquid shift of water.

Conclusion

In conclusion, we present in this work a thorough analysis of the evaluation of the 17 O and 1 H gas-to-liquid shifts of water at the DFT level. We consider all the relevant parameters that can impact such property: basis set, xc-functional, structural model, statistical convergence and approach to calculate nuclear shielding. Although both molecular and periodic approaches provide important insight into the correct way to calculate gas-to-liquid shifts, we demonstrate that both approaches does not behave in the same way when basis set and xc-functional are considered.

From the periodic calculations, which can be considered as converged in terms of statistical sampling and basis-set extension and that do not suffer from finite-size effects, the present results show that the 17 O gas-to-liquid shift is extremely sensitive to the xc-functional used for the GIPAW evaluation of the nuclear shieldings whereas the impact of the structural model is moderate. We have shown that it results from strong differences in the description of electronic density de-localization between water molecules: PBE and other GGA xc-functionals tend to over-delocalize the electronic density from the oxygen atoms to the adjacent hydro-gen atoms while the inclusion of HF exchange correct this shortcoming. Nevertheless, an xc-functional including too much HF exchange lead to a too strong re-localization of the density on the oxygen atoms which is as much prejudicial for gas-to-liquid shift estimates.

Comparisons between six different xc-functionals show that PBE0 and B3LYP provide a good compromise to correctly describe the electronic densities. It is worth mentioning that such dependency of the 17 O gas-to-liquid shift towards local electronic de-localization prevents the use of an implicit description of the solvation shell around a given water molecule.

This point is further demonstrated and discussed in the Supporting Information. Due to a narrower nuclear shielding range, the impact of the xc-functional on 1 H gas-to-liquid shift is moderate, 1.6 ppm at most with the six considered xc-functionals. However, none of them leads to the experimental value which suggests that there is still some room for improvements, starting from the present simulations. We can assume the NQEs are likely to have a significant impact on this value.

Molecular calculations are more difficult to rationalize. Indeed, despite careful analysis of the error originating from finite-size effects and statistical sampling, contributions originating from basis-set effects, xc-functional and structural model remain intertwined. As for GIPAW calculations, the description of the electronic density de-localization between water molecules is important. Fortuitously, the use of a small basis set limit the de-localization of the density which may lead to a 17 O gas-to-liquid close to the experimental value. In contrast, extended basis set are enable to provide values close to the experimental gas-to-liquid shift. The 1 H gas-to-liquid appears less influenced by basis set effects.

The present study shows that the [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H gas-to-liquid shifts of water are properties which accurate calculation is none trivial. However, we provide a detailed guideline to properly evaluate them. Such guideline can be extended to other hydrogen-bonded liquids such as ammonia, or even for solvated species in highly polar solvents.

The present study shows that the [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H gas-to-liquid shifts of water are properties which accurate calculation is none trivial. This is due to the strong hydrogen bonds that exist between the water molecules and that lead to electronic density delocalization which is difficult to describe theoretically. However, we provide a detailed guideline to properly evaluate those two properties. Such guideline can be extended to other hydrogen-bonded systems such as liquid ammonia and the evaluation of its 14 N gas-to-liquid shift. In this latter case, the existence of weaker hydrogen bonds between NH 3 molecules could suggest more limited basis-set effects in molecular calculations and a weaker impact of the xc-functional.

The methodology could also be extended to species solvated in water, such as H 3 PO 4 which is the experimental reference for 31 P NMR spectroscopy.

Supporting Information Available

Discussion on the computational procedure applied to calculate nuclear shieldings in the gas phase and accuracy of implicit solvent model. Additional test on basis-set convergence of NMR GIAO calculations. Distribution of [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O and 1 H average nuclear shielding per atom obtained over the 50 considered MD frames. Gas-to-liquid shifts obtained from the MD-PBE-300 molecular dynamics simulation. GIPAW gas-phase reference nuclear shielding values used to calculate the gas-to-liquid shifts presented in the main document.

Figure 1 :

 1 Figure 1: (Left panels) 17 O (top) and 1 H (bottom) nuclear shieldings convergence with respect to the cut-off energy in GIPAW calculations for a single water molecule in the gas phase at the PBE (black circles) and PBE0 (red circles) levels. (Right panels) Similar data for one given molecule (circles) in a randomly chosen MD frame as well as for the mean value over the 64 water molecules (square) of the same frame. Both PBE (black) and PBE0 (red) xc-functionals were also used. The dashed lines are guide for the eyes to indicate fully converged values.

Figure 2 :

 2 Figure 2: 1 H (top) and 17 O (bottom) nuclear shielding convergence of the central water molecule in H 2 O-(H 2 O) n clusters with respect to n in NMR GIAO calculations. Each value (red dot) was obtained from an average over 1025 configurations (25 distinct central water molecules along 41 distinct frames) extracted from the MD-PBE-300 trajectory. The plain black bars represent the standard deviations. All calculations were performed at the B3LYP/6-311G(d,p) level.

and 1 .

 1 56 ppm for 1 H and[START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O, respectively, and 0.23 and 2.11 ppm when only half the frames are considered. Those values are rather low considering the respective nuclear shielding ranges and can be considered as a good measure of the statistical error that exist in the results presented below. Although all the aforementioned calculations were performed at the B3LYP/6-311G(d,p) level of theory only, similar conclusions on the statistical convergence would be obtained with other combinations of xc-functional and basis set.In the light of those preliminary results, in the section , we generated 1000 water aggregates, obtained from 50 different frames randomly extracted from MD calculations, each one centered on one water molecule, surrounded by around 20 (± 3) other water molecules that are located up to 6 Å from it. Only nuclear shieldings of the central molecule were extracted in order to compute, using the NMR GIAO approach, mean values and standard deviations presented and discussed in Section .

Figure 3 :

 3 Figure 3: 1 H (top) and 17 O (bottom) nuclear shielding convergence of the central water molecule in H 2 O-(H 2O) 36 as a function of the number of configurations taken into consideration in NMR GIAO calculations. All configurations were extracted from the MD-PBE-300 trajectory. Three quantities are considered: means (black lines), standard deviations (red lines) and standard deviations over the mean (blue lines) for an increasing number of configurations. All calculations were performed at the B3LYP/6-311G(d,p) level.

  Four xc-functionals, namely PBE, PBE0, B3LYP and HSE, in combination with eight different basis sets of increasing size were applied to calculate 1 H and 17 O gas-to-liquid shifts.The results are presented in Figure4. Various observations can be done. First, regardless of the xc-functional, 6-31G(d,p) and 6-31G(2d,2p) systematically overestimate the 17 O gasto-liquid shift while the six other basis sets always underestimate it. Second,[START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O gas-toliquid shift increases from the smallest basis sets to the largest ones, from ∼-25 ppm for 6-31G(d,p) to ∼-57 ppm for 6-311++G(2d,2p). This behavior is the same whatever the considered xc-functional. Furthermore, Figure4also demonstrates that diffusion functions on the oxygen atoms are crucial in the calculation of the 17 O gas-to-liquid shift. Indeed, a difference of almost 25 ppm is obtained between 6-31G(d,p) and 6-31+G(d,p) results. In contrast, adding polarisation functions or going from double-to triple-ζ has a more limited impact. Interestingly, adding diffusion functions on the hydrogen atoms has only a minor influence on the 17 O gas-to-liquid shift. Regardless of the basis set, hybrid xc-functionals (B3LYP, HSE and PBE0) lead to larger gas-to-liquid shifts as compared to PBE. From 6-

Figure 4 :

 4 Figure 4: 17 O (Top) and 1 H (Bottom) theoretical gas-to-liquid shift as a function of the basis set used for NMR GIAO calculations. Four xc-functionals have been considered: PBE (black), B3LYP (red), HSE (blue) and PBE0 (green). The black dashed line represents the experimental value. 9-13 All calculations were performed on a single H 2 O-(H 2 O) 42 water cluster extracted from the MD-PBE-300 simulation. Additional data are presented in Figure S1 to S3 of the Supporting Information.

Figure 5 :

 5 Figure5: The oxygen-oxygen (g OO (r), bottom) and oxygen-hydrogen (g OH (r), top) radial distribution functions of liquid water obtained from the MD-PBE-300 (black), MD-OptB88-300 (green) and MD-PBE-330 (red) molecular dynamics simulations. The experimental data of Soper are also provided for comparison.[START_REF] Soper | The Radial Distribution Functions of Water as Derived from Radiation Total Scattering Experiments: Is There Anything We Can Say for Sure[END_REF]

O and 1

 1 H gas-toliquid shifts calculated by means of molecular approach for the three considered MD simulations: MD-PBE-300, MD-PBE-330 and MD-OptB88-300 are presented in Figure6. Four basis sets were used: 6-31G(d,p), 6-31G(2d,2p), 6-311G(d,p) and 6-311++G(d,p) for the MD-PBE-300 trajectory and only 6-31G(d,p) and 6-31G(2d,2p) for the MD-PBE-330 and MD-OptB88-300 trajectories. Considering MD-PBE-300 only, 6-31G(d,p) and6-31G(2d,2p)

Figure 6 :

 6 Figure 6: 17 O (Top) and 1 H (Bottom) difference between experimental and theoretical gas-toliquid shift for three distinct molecular dynamics simulations (MD-PBE-300, MD-OptB88-300 and MD-PBE-330). All GIAO calculations are performed using the B3LYP xc-functional. Four different basis sets are employed: 6-31G(d,p) (blue bar), 6-31G(2d,2p) (red bar), 6-311G(d,p) (black bar) and 6-311++G(d,p) (green bar). Error bars are equal to ±0.23 and ±2.11 ppm for 1 H and 17 O, respectively, and are evaluated from the standard deviations over the mean discussed in section .

Figure 7 Figure 7 :

 77 Figure7displays the distributions of 1 H and[START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shielding values obtained at the PBE and PBE0 levels using 50 frames extracted from the MD-OptB88-300 trajectory. Each distribution is well defined which further confirms that a correct statistical sampling of the structural features is achieved. In FigureS4of the Supporting Information, we also present the distribution of the average nuclear shielding per atom obtained over the 50 considered frames. The four distributions are much narrower than the distributions of Figure7. Indeed, standard deviations extracted from Figure7are 2.2 (1.9) and 12.8 (8.0) ppm at the PBE (PBE0) level for 1 H and 17 O, respectively, while they are equal to 0.9 (0.7) and 4.8(3.1) ppm in FigureS4of the Supporting Information. This demonstrates that the spread of the distributions presented in Figure7are due to thermal fluctuations of the local environment and not to a distribution of chemical environments undergone by each water molecule.

Figure 8 :

 8 Figure 8: (Top) Gas-phase H 2 O isosurface decomposition of the electronic density difference (∆ρ) between PBE0 and PBE calculations. (Bottom) Same decomposition applied to a single MD snapshot of liquid water. Colors are used to indicate the sign of ∆ρ: red for positive values and blue for negative ones. The red isosurface means +0.001 e/ Å3 , when the blue is -0.001 e/ Å3 .

  [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O nuclear shielding range. Consequently, we can consider that clusters containing at least 16 H 2 O molecules leads to converged nuclear shielding values for both hydrogen and oxygen.

			32
			31
			30
	σ H iso [ppm]	26 27 28 29
			25
			24
			0 23	4	8 12 16 20 24 28 32 36
		340
		330
		320
	[ppm]	310
	σ O iso	300
		290
		280
		0 270	4	8 12 16 20 24 28 32 36
				Cluster Size (n)

  [START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF] O are obtained between MD-PBE-300 and MD-OptB88-300. Based on the MD-PBE-330 simulation in combination with the 6-31G(2d,2p) basis set, the experimental 17 O gas-to-liquid shift is almost recovered: -33.8 vs -36.1 ppm. Similarly, the best GIAO estimate for ∆σ iso H equals -4.7 ppm and is obtained for the same MD simulation in combination with the 6-31G(d,p) basis set. In both cases,

Table 1 :

 1 Mean gas-to-liquid shifts obtained from 50 randomly chosen snapshots taken from the MD-OptB88-330 and MD-PBE-330 molecular dynamics simulations. Bold notations (PBE, optB88-vdW, HSE, B3LYP, PBE0 and PBE0-50%) denote the xc-functional used for the GIPAW calculations. All values are given in ppm, when standard deviations are given in parenthesis.

			MD-OptB88-300			MD-PBE-330		Exp.
		PBE	optB88-vdW	HSE	PBE	optB88-vdW	HSE
	1 H	-7.0 (2.2)	-7.0 (2.3)	-6.1 (2.0)	-7.2 (2.6)	-	-6.2 (2.2)	-4.3 a
	17 O -64.7 (12.8)	-64.2 (12.8)	-46.6 (8.6)	-65.8 (14.4)	-	-47.5 (9.7)	-36.1 b
		B3LYP	PBE0	PBE0-50%	B3LYP	PBE0	PBE0-50%
	1 H	-6.3 (2.0)	-6.1 (1.9)	-5.4 (1.5)	-6.3 (2.2)	-6.2 (2.0)	-5.4 (1.8)	-4.3 a
	17 O -42.3 (8.8)	-37.6 (8.0)	-23.1 (5.5)	-43.1 (10.0)	-38.3 (9.1)	-23.6 (6.2)	-36.1 b
	a From Ref. 9,12
	b From Ref. 10,11,13

Table 1 ,

 1 one can see that the six considered xc-functionals are not equivalent to describe such interaction, especially when looking at ∆σ iso O values. The two GGA xc-functionals, PBE and optB88-vdW, lead to the worst results with an error of almost 100% as compared to experiments. PBE0-50% is the only xc-functional leading to a value that is smaller in absolute values than

Acknowledgement

This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2017-p1320. This work was also granted access to the HPC resources of CINES under the allocation 2017-A0020810077 made by GENCI. The authors declare that there has been no significant financial support for this work.