{'fraction_non_alphanumeric': 0.025858093418538772, 'fraction_numerical': 0.019830028328611898, 'mean_word_length': 4.787075867327488, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Static and dynamical properties of heavy water at ambient conditions from first-principles molecular dynamics

The static and dynamical properties of heavy water have been studied at ambient conditions with extensive Car-Parrinello molecular-dynamics simulations in the canonical ensemble, with temperatures ranging between 325 K and 400 K. Density-functional theory, paired with a modern exchange-correlation functional (PBE), provides an excellent agreement for the structural properties and binding energy of the water monomer and dimer. On the other hand, the structural and dynamical properties of the bulk liquid show a clear enhancement of the local structure compared to experimental results; a distinctive transition to liquid-like diffusion occurs in the simulations only at the elevated temperature of 400 K. Extensive runs of up to 50 picoseconds are needed to obtain well-converged thermal averages; the use of ultrasoft or norm-conserving pseudopotentials and the larger plane-wave sets associated with the latter choice had, as expected, only negligible effects on the final result. Finite-size effects in the liquid state are found to be mostly negligible for systems as small as 32 molecules per unit cell.


I. INTRODUCTION
Water, due to its abundance on the planet and its role in many of the organic and inorganic chemical processes, has been studied extensively and for decades both at the theoretical and at the experimental level 1−17 . The peculiar interplay of hydrogen bonding, glassy behavior, and of quantum-mechanical effects on the dynamics of the atomic nuclei make computer simulations challenging, and a great effort has been expended to build a comprehensive and consistent microscopic picture, and a link with observed macroscopic properties 4−17 . Additionally, it is only recently that close agreement for fundamental structural information such as the radial distribution function has been obtained between different experimental techniques, such as X-ray 2 and neutron diffraction 3 measurements.
Computational studies based on molecular dynamics simulations have also a rich history in the field. Simulations using force-fields models 4−9 have been successful at reproducing many structural and dynamical properties of liquid water. However, empirical models rely on parameters which are determined by fits to known experimental data, or occasionally to ab-initio results. Their transferability to different environments, or the ability to reproduce faithfully the microscopic characteristics of hydrogen bonding, are often in question. Due to development of novel techniques 18,19 and the ever-increasing improvement in computational power, extensive moleculardynamics simulations from first-principles are now possible. The increased accuracy and predictive power of these simulations comes at a significant price, and careful considerations has to be given to the length scales and time scales that can be afforded in a first-principle simulation, and the trade-offs in statistical errors when compared with classical simulations. Numerous ab-initio simulations on water have appeared 10−17 , showing good agreement with experiments for the structural and dynamical data. Recent results have also reported 20,21 that careful equilibration for ab-initio water at ambient conditions leads to radial distribution functions over-structured in comparison with experiments 2 . After equilibration, the numerical estimates for the diffusion coefficient become at least one order of magnitude smaller than the measured ones.
Prompted by these results and by our own observations on the structure of water around iron aqua ions 22 , we have undertaken an extensive investigation of the static and dynamical properties of water, to ascertain its phase stability around ambient conditions as predicted by firstprinciples molecular dynamics. Particular care has been given to the statistical accuracy of the results, assuring that the time scales and length scales of the simulations were chosen appropriately for the given conditions. The paper is organized as follows: In section II, we detail all the technical aspects of our simulations. Section III surveys the static and vibrational properties of the water molecule and the water dimer in vacuum, at the GGA-PBE 23 density-functional level. In section IV, we discuss our extensive liquid water simulations, performed with Car-Parrinello molecular dynamics, in the temperature range between 325 K and 400 K. Section V discusses the limitations of this approach, and some of the possible reasons to explain the remaining discrepancies with experimental results.

II. CAR-PARRINELLO MOLECULAR DYNAMICS
Our first-principles calculations are based on densityfunctional theory, periodic-boundary conditions, planewave basis sets, and norm-conserving 24 or ultrasoft pseudopotentials 25 to represent the ion-electron interactions, as implemented in the public domain codes CP and PWSCF in the ν-ESPRESSO package 26 .
In Car-Parrinello molecular dynamics, an extended Lagrangian is introduced to include explicitly the wavefunction degrees of freedom, that are evolved "on-the-fly" simultaneously with the ionic degrees of freedom: In Eq. 1, Ψ i (r) are the occupied Kohn-Sham orbitals, µ is a fictitious mass parameter used to control the evolution of the electronic degrees of freedom in time, M I are ion masses, E KS is the Kohn-Sham energy and Λ ij are Lagrange multipliers, used to impose the orthonormality constraint Ψ * i Ψ j = δ ij . An extensive review of the method can be found in Ref. 27 ; the subtle technical issues arising in the simulations of liquid water have been discussed in the recent literature 20,28 . We note in passing that our simulations' parameters agree with the recommendations set forth in these latter papers. It should be stressed that full convergence of the forces acting on the nuclei with respect to the basis set is easily achieved with plane waves, and no Pulay forces arise in the dynamics, since the basis set is independent from the atomic positions.

A. Technical Details
The structural and vibrational properties of the water monomer and dimer and the binding energy of the dimer have been calculated using density-functional theory in the generalized-gradient approximation and the total energy pseudopotential method, and densityfunctional perturbation theory 29 , as implemented in PWSCF. We performed separate calculations using either norm-conserving pseudopotentials for both the hydrogen and the oxygen, or ultrasoft ones. These same pseudopotentials were also used for the norm-conserving or ultrasoft molecular dynamics simulations. In particular, the O Troullier-Martins norm-conserving pseudopotential was generated using the FHI98PP package 30 with core radii for the s, p and d components of 1.25 a.u., 1.25 a.u., and 1.4 a.u. respectively. The Troullier-Martins hydrogen pseudopotential was generated using the Atom code 31 with a core radius for the s component of 0.8 a.u. The ultrasoft pseudopotentials were taken from the standard PWSCF distribution 32 . The Kohn-Sham orbitals and charge density have been expanded in plane waves up to a kinetic energy cutoff of 25 and 200 Ry (respectively) for the ultrasoft case, and of 80 Ry and 320 Ry for the norm-conserving case. A cubic supercell of side 30 a.u. was used; interaction with periodic images is negligible 33 with this unit cell size.

III. WATER MONOMER AND DIMER: STRUCTURAL AND VIBRATIONAL PROPERTIES
The equilibrium structures and energetics are summarized in Table I. We have included published results 33 using the BLYP functional for comparison. Both ultrasoft and norm-conserving PBE density functionals show very good agreement with experimental values. In particular, the PBE results have a dimer binding energy in closer agreement to the experiments than BLYP; the binding energy in this latter case is too weak by 4 kJ/mol, and exhibits a longer O-O distance. Table II and III show respectively the vibrational frequencies of the water monomer and dimer in vacuum. In this calculation, a hydrogen mass of 1 a.m.u. was used (as opposed to the 2 a.m.u. mass used for the dynamical simulations of heavy water). The calculations of the vibrational frequencies with PWSCF using a cubic cell of size (30 a.u.) 3 . To achieve a convergence of a few cm −1 in the frequencies, cutoffs of 35 Ryd and 420 Ryd were used for wavefunctions and charge densities, respectively, in the ultrasoft case and, in the norm-conserving case, 100 Ryd and 400 Ryd were used. As shown in Table II and III,the PBE functional gives intramolecular stretching modes that are in general blue-shifted compared to BLYP and experimental results. In the calculations for the dimer, the libration modes are also higher for the PBE functionals then those given by experiments and BLYP. We note in passing that the errors on these frequencies (especially the low energy ones) are slightly larger than usually expected from DFT). We will return  to this point in later section.

IV. LIQUID WATER SIMULATIONS
A. Liquid water simulation at 325 K

Simulation Details
In this first simulation, we used a body-centered-cubic supercell with 32 heavy water molecules, periodic boundary conditions, and the volume corresponding to the experimental 40 density of 1.0957 g/cm 3 at 325 K. A body-centered-cubic supercell strikes the optimal balance, for a given volume, in the distance between a molecule and its periodic neighbors, and the number of these periodic neighbors. Ultrasoft pseudopotentials were first used, as detailed in the previous section, with plane-wave kinetic energy cutoffs of 25 Ry (wavefunctions) and 200 Ry (charge densities). The deuterium mass was used in place of hydrogen to allow for a larger time step of integration. It should be noted that for classical ions this choice does not affect thermodynamic properties such as the melting temperature (the momentum integrals for the kinetic energy factor out in the Boltzmann averages). Of course, dynamical properties such as the diffusion coefficient will be affected by our choice of heavier ions. Extensive experimental data for deuterated (heavy) water are in any case widely available.
The wavefunction fictitious mass (µ) is chosen to be 700 a.u.; this results in a factor of ∼14 between the average kinetic energy of the ions and that of the electrons. A time step of 10 a.u. was used to integrate the electron and ionic equations of motions. This combined choice of parameters allows for roughly 25 ps of simulation time without a significant drift in the kinetic energy of wavefunctions and the constant of motion for the Lagrangian (1). Our choice of fictitious mass is consistent with the ratio µ/M ≤ 1 3 for heavy water molecules suggested by Grossman et al. 20 , and assures that the physical properties are not influenced by the electronic degrees of freedom. Our initial configuration was obtained from a comparatively short 1.2 ps simulation at twice the value of the target temperature (650 K instead of 325 K); a restart with zero initial velocities was then performed at 325 K, with the temperature controlled by a single Nose-Hoover thermostat on the ions (no electronic thermostat is applied in any of the simulations).

Results
The thermostat stabilizes quickly (∼1 ps) the temperature around the target value of 325 K, with the usual fluctuations due to small size of the system. However, the system is still far from equilibrium; this can be clearly observed by looking at the time evolution of the radial distribution function (RDF) and mean square displacements (MSDs). We plot in Fig. 1 the MSDs of the oxygen atoms as a function of time. For the first 10 ps, the water molecules diffuse with a velocity comparable with experimental data; after 10 ps, a sharp drop in diffusivity is observed, accompanied by a distinctive sharpening of the features in the oxygen-oxygen radial distribution function (see Fig. 2). The radial distribution functions was calculated from the infinite bulk system by repeating the unit cell in all directions. This means that molecules up to 5.5Å in a 32-molecule cell are inequivalent. As we move beyond this cutoff distance, the radial distribution functions will include both molecules that are inequivalent, and some that are equivalent. The statistical accuracy is going to be gradually, but slowly, affected. In this graph, and most of the latter figures, we plot radial distribution functions only up to 5.5Å(but see e.g. Fig. 10 for a discussion of the finite-size effects).
The potential energy for the dynamics (Fig. 1) also drifts downward in these first 10 ps, and stabilizes afterward. We calculated the self-diffusion coefficient (D self ) from the Einstein relation (in 3 dimensions): The structural and dynamical properties before and after this 10 ps mark are summarized in Table IV; experimental values at 298 K are included for comparison. All these observation conjure to a picture in which the system takes at least 10 ps to reach a reasonably thermalized  state, in a process somewhat reminiscent of a glass transition. Although the time needed for equilibration will be dependent on the initial conditions, these preliminary result suggests that simulation times in the order of ten of picoseconds might be need to calculate well-converged thermodynamic observables.
Once the initial thermalization trajectory was discarded from our averages, we obtained a self-diffusion coefficient one order of magnitude smaller than the experimental value measured at room temperature. This result, combined with the clear over-structuring of the oxygen-oxygen radial distribution function g OO (r), indicates that our system has reached a "frozen" equilibrium state very different from what expected for liquid water (strictly speaking, a system with a finite and small number of inequivalent atoms or molecules will never undergo a phase transition).
These considerations indicates the need to assess accurately the phase stability of liquid water as obtained from first-principle molecular dynamics. At the same time, they point out the requirement of long simulation times, and a careful analysis of the technical details of the simulations. B. Extensive water simulations in the region between 325 K and 400 K

Simulation details
In order to find out the temperature at which our system would move from a glassy to a liquid-like state, we decided to perform a series of extensive simulations at increasing temperature, from 325 K to 400 K. We first performed ≈ 25 ps simulations at 325 K, 350 K, 375 K and 400 K, using µ=700 a.u. and δt=10 a.u. . We used at every temperature the experimental densities 40 , with the caveat that the 400 K value was obtained by extrapo- lation. 25 ps is approximately the maximum time allowed for a simulation with these parameters before the drift in the kinetic energy of the wavefunctions becomes apparent. We thus used these simulations as efficient "thermalization" runs, to be followed by production runs that will be described below. It is interesting to monitor during these thermalization runs the evolution of the MSDs; these are shown for all four temperatures in Fig. 3. An abrupt drop in diffusivity is observed for all cases but one (predictably, the one at the highest temperature of 400 K). The onset of this drop in diffusivity varies, but broadly speaking is again of the order of 10 ps. With these trajectories, now well thermalized in configuration space, we started our four production runs at 325 K, 350 K, 375 K and 400 K, and each of them starting from the last ionic configurations of the previous simulation at the corresponding temperature, but with zero ionic velocities. These four runs, lasting between 20 ps and 37 ps, were performed using µ=450 a.u. and δt=7 a.u. This choice of mass and timestep allows for an excellent conservation of the constant of motion, and negligible drift in the fictitious kinetic energy of the electrons, for the simulation times considered. The ratio between the kinetic energy of the ions and that of the electrons was ≈ 22 for the whole production time. Although a small enough fictitious mass decouples the electronic and ionic degrees of freedom, Tangney et. al. 41 pointed out that there is a fictitious mass dependent error that is not averaged in the time scale of ionic motions. Schwegler et. al. 28 studied this effect comparing closely Car-Parrinello and Born-Oppenheimer simulations finding a larger self-diffusion coefficient in the Car-Parrinello simulation. However, the structural and thermodynamical properties were not affected. We show in Figs. 4 and 5 the case of the 400 K simulation; we stress that no periodic quenching of the electrons was needed, and the simulations were single uninterrupted runs. Since the initial configurations were already at equilibrium at their respective temperature, and thermalization in momentum space is fast, we found that "production time" can start early in the simulations. We discarded from each trajectory the initial 1.2 ps that were needed to allow the ions to reach their target kinetic energy. As a measure of the good thermalization reached in the simulations, we show in Fig. 6 the O-O radial distribution function obtained from the first 12 ps of our 400 K trajectory, and the following 12 ps.
To rule out any spurious effect in our simulations coming from the use of pseudopotentials, or an extended La- grangian, we also performed two simulations using normconserving pseudopotentials (as described in Section II). These require larger plane-wave basis sets (80 Ry for the wavefunctions and 320 Ry charge densities, corresponding to 40000 plane waves vs 7000 for the ultrasoft case), and discrepancies, if any, with the ultrasoft calculation will provide an approximate estimate of the effects of the pseudopotential approximation and of the dynamics of the fictitious degrees of freedom. For the runs involving these norm-conserving pseudopotentials, we used µ =300 a.u. and δt=5 a.u 42 . These parameters results in a factor of ∼ 13 between the kinetic energy of the ions and that of the electrons. Details of all these simulations are shown in Table V.

Results
The oxygen-oxygen and oxygen-deuterium radial distribution functions for the different conditions considered in this work are shown in Fig. 7 and Fig. 8 respectively. We find that at temperatures of 375 K and below both g OO and g OH show considerably more structure than found experimentally at 300 K. For this range of temperatures the height of the first peak of g OO (r) is roughly between 3.2 and 3.4, and significantly larger than the experimental value of 2.75 (also measured at 300 K). However, when the temperature in our simulations is increased to 400 K, a distinct drop of the first peak to ∼ 2.5 is observed, the radial distribution function g OO (r) and g OH (r) show a sharp change in their structure, and the water molecules start diffusing much faster, as reflected in the MSDs curves for the oxygen atoms shown in Fig. 9. To provide cleaner statistics, the MSDs curves shown have been calculated as an average over individual MSDs curves, each obtained from our trajectory by shifting -for each individual MSDs curve -the starting configuration by 0.017 ps (in this way, a 17 ps trajectory would provide 1000 progressively shorter MSDs curves that are then averaged). The self-diffusion coefficient D self is calculated from the slope of the respective MSDs curve in the range of 1 to 20 ps using Einstein's relation (2). Negligible differences are observed between simulations performed with ultrasoft or norm-conserving pseudopotentials, ruling out any role of the pseudopotential details in this observed result. The structural and dynamical results are summarized in Table VI. As seen in this table, there is an eight-fold increase in D self when increasing the temperature from 375 K to 400 K. Price et. al. 43 reported the experimental self-diffusion coefficient of supercooled heavy water at different of temperatures. At 276.4 K, which is just below the freezing temperature of heavy water (277.0 K), the experimental value for D self is 0.902 × 10 −5 cm 2 /s. D self for our simulations at 325 K, 350 K and 375 K is 0.16, 0.25 and 0.26 ×10 −5 cm 2 /s, respectively. These numbers are significantly smaller then the experimental value below the freezing point. On the other hand, D self at 400 K in our simulations is comparable to the experimental value at 300 K. These observations suggest that the theoretical freezing point for water at the DFT-PBE level is between 375 K and 400 K, and water below 375 K is in a glassy/supercooled state.
The hydrogen-bond structure can be studied calculating the number of hydrogen bonds per molecule: we identify a hydrogen bond is identified when two oxygen atoms are closer than 3.5Å and the OHO angle is greater than 140 o (consistently with Ref 44 , and at slight variance with Ref 45 ). The results are shown in the last column of Table VI. Between the temperatures of 325 K and 375 K there are only small changes in the number of hydrogen bonds in the system, going from 3.86 per molecule to 3.72 per molecule. An abrupt decrease to 3.37-3.45 bonds occurs when the temperature increases from 375 K to 400 K. The experimental value of 3.58 at 300 K lies between our values at 375 K and 400 K, further suggesting that the freezing point in our simulations is between 375 K and 400 K.
In summary, we found clear liquid-like signatures in the structure and dynamics of (heavy) water, as described by DFT-PBE and Car-Parrinello MD, for temperatures reaching at least 400 K. At temperatures of 375 K and lower water is found to be in a glassy state, more structured and with much lower diffusivity. This discrepancy of more than 100 K between experimental and theoretical results is obviously relevant, given the enormous importance of water in the description of systems ranging from electrochemistry to biology, and it is investigated further in the next section.

V. OVERESTIMATION OF THE FREEZING TEMPERATURE OF WATER
There are several possible reasons for the overestimation of the freezing temperature of water, and several of them could play a significant role. We discuss here some of the possibilities.

A. Finite-size effects
Since our simulation cell contains only 32 water molecules, finite-size effects could obviously play a role even if periodic-boundary conditions are used. The interactions of water molecules with their periodic images could be considerable due to the long-range hydrogenbond network. On the other hand, when zero correlations are found between a molecule and its periodic image we can safely assume that the unit cell is for all practical purposes large enough, and every molecule feels the same environment that it would have in an infinite system. In our case, the distance between a molecule and its eight periodic images is ∼ 11Å and at this distance all radial distribution functions look very flat and unstructured. In any case, to study the finite-size effects we carried out another extensive simulation (40 ps total, with 15 ps of production time following 25 ps of thermalization) for a system composed of 64 heavy water molecules at 400 K. We used the same parameters for this simulation as in the 32-molecule, 400 K ultrasoft simulation. The oxygen-oxygen radial distribution function g OO (r) is shown in Fig. 10 for both the 32-and 64-molecule systems. As mentioned previously, the radial distribution functions was calculated by repeating the unit cell in all directions. The molecules up to 5.5Å in the 32-molecule cell, and 6.9Å in the 64-molecule cell are inequivalent. We indicate in the graphs with two arrows, the radii of the spheres completely inscribed by in BCC simulation cells.
Differences between 32-and 64-molecule systems are negligible, and within the variance for simulations of the order of 10-20 ps (as estimated from uncorrelated classical simulation data 20 ); the 64-water simulation shows  a marginally more structured g(r) where the first peak height is at about 2.6 (compared to 2.5 for the 32-water case). Larger ab-initio simulations would be too demanding; for this reason, we performed two classical simulations at 300 K using the SPC force field 47 for water, and comparing the case of 64 and 1000 water molecules (each simulation lasting 1000 ps). The g OO (r) calculated from these two runs are shown in Fig. 11, and again we do not find any significant differences between these two curves. These results help ruling out finite-size effects as the major cause of the discrepancy observed with the experimental numbers.

B. Exchange-correlation functional effects
While density-functional theory is in principle exact, any practical application requires an approximated guess to the true exchange-correlation functional. In this work, we have used the GGA-PBE approach 23 . As it was observed in Sec. III, the structural properties for the water molecule and dimer are in excellent agreement with experiments, as is the binding energy for the dimer. On the other hand, the vibrational properties show larger discrepancies with experiments than usually expected, in particular for some of the libration modes in the dimer. This result certainly points to the need for improved functionals to describe hydrogen bonding. The dependence of the melting point on the exchange-correlation functional chosen is more subtle; below 400 K, PBE water displays solid-like oxygen-oxygen radial distribution functions that are only slightly affected by the temperature, and that are similar to those obtained with a fairly different functional such as BLYP 20 . The similarity between these radial distribution functions is just a reminder that the structural property and the geometry of the intermolecular bonds are well described by different functionals; once water is "frozen", all radial distribution functions will look similar. The temperature at which this transition takes place could be affected by the use of different functionals 48 and the magnitude of the contribution of any one of them to the melting point temperature is still open to investigation.

C. Quantum effects
In first-principles Car-Parrinello or Born-Oppenheimer molecular dynamics simulations the ions are most often treated as classical particles, which is a good approximation for heavy ions (path-integral simulations can describe the quantum nature of the ions 49−51 , but their computational costs, when paired with a firstprinciples DFT descriptions of the electrons, preclude at this moment simulations with the statistical accuracy needed). However, for light ions like hydrogen or deuterium, the effects of a proper quantum statistics can be very significant; tunneling of the nuclei can also affect the dynamics 52,53 . In the case of water, all the intramolecular vibrational modes and some of the intermolecular modes are much higher in energy compared to room temperature. We show in Fig. 12 the power spectrum for the deuterium atoms as calculated from the velocity-velocity correlation function of heavy water molecules for our simulations at 400 K (ultrasoft, 32 molecules). The distinctive peaks of the intramolecular stretching and bending modes are centered around 2400 cm −1 and 1200 cm −1 respectively, much higher than the room temperature of k B T room ≈ 200 cm −1 . The peak at 500 cm −1 corresponds to the intermolecular vibrational modes, also larger than k B T room . When ions are treated as classical particles, as in our ab-initio molecular dynamics simulations, all vibrational modes obey Boltzmann statistics. In reality, modes with frequency higher than k B T room are frozen in their zero-point motion state, and their exchange of energy with the lower-frequency modes ("the environment") is suppressed -in other words their contribution to the specific heat is zero, in full analogy with the lowtemperature discrepancies from the Dulong-Petit law in solids. This effect could significantly affect the dynamics of water-water interactions, and it has long been argued that treating each water molecule as a rigid body could actually provide a closer match with experimental conditions. In fact, recent ab-initio simulations 54 in which the water molecules are constrained to maintain their equilibrium intramolecular bond lengths and bond angle result in a more diffusive and less structured description of water, that remains liquid at a temperature of 326 K. Path-integral simulation 55 for water described with classical potentials did find as well a significant difference due to the quantum effects (i.e. the freezing of the highenergy vibrational excitation), of the order of 50K.

VI. CONCLUSION
We performed extensive first-principles molecular dynamics simulations of heavy water at the DFT-PBE level. Equilibration times are found to be comparatively long, and easily in excess of ≈ 10 ps at ambient temperature. Good statistics was obtained with simulations that included at least 25 ps of thermalization time, followed by 20 to 40 ps of production time. At ambient temperature, water is found to be over-structured compared to experiment, and with a diffusivity that is one order of magnitude smaller than expected. An abrupt change in the structure and dynamics is observed when the the temperature is raised from 375 K to 400 K; even at this high temperature, where a liquid-like state is reached, water shows more structure and less diffusivity than found experimentally at room temperature. These results are broadly independent of some of the possible errors or inaccuracies involved in first-principles simulations, including in this case finite-size effects, insufficient thermalization or simulation times, or poorly designed pseudopotentials. Our simulations suggest that the freezing point is around 400 K -this discrepancy is at variance with the good agreement for the structure and energetics of the water monomer and dimer with the experimental values, and could originate in the neglect of quantum statistics for the many high-frequency vibrational modes in this system. This work was performed under the support of the Croucher Foundation and Muri Grant DAAD 19-03-1-0169. Calculations in this work have been done using the ν-ESPRESSO package 26