Ensemble first-principles molecular dynamics simulations of water using the SCAN meta-GGA density functional

Molecular dynamics (MD) simulations of water have been the subject of numerous studies for nearly 50 years, beginning in the 1970s with the work of Rahman and Stillinger.^1–3 Modeling of water using classical molecular dynamics has been refined over time with the development of increasingly accurate force fields, an area of research that is still very active.^4–9 First-principles quantum mechanical simulations have become more readily available with the degree of computing power available today and with the development of density functional theory (DFT)^10–12 coupled with ab initio molecular dynamics¹³ (MD) simulations.

The accuracy of DFT simulations has progressively improved over time as new exchange-correlation functionals have been introduced.¹⁴ The simplest density functional is the local density approximation (LDA) functional¹⁵ that does not accurately describe hydrogen bonds, leading to poor chemical accuracy. Accuracy can be somewhat improved by using the generalized gradient approximation (GGA).¹⁶ In particular, several authors have found that the GGA functional developed by Perdew, Burke, and Ernzerhof (PBE)¹⁷ can provide a reasonable description of liquid water,^18,19 although PBE simulations need to be run at an elevated temperature (e.g., 400 K) in order to approximate the properties of the liquid at 300 K. A large number of more elaborate density functionals have been proposed, notably hybrid density functionals that include a fraction of the Hartree-Fock exchange energy,^20–22 van der Waals functionals^23–26 that include a description of long-range dispersion interactions, and meta-GGA functionals^27–29 that involve a kinetic energy density term. When used in simulations of water, the addition of a van der Waals term to a simpler functional such as PBE has been shown to improve the chemical accuracy of the functional by correcting the intermediate range intermolecular interactions in liquid water leading to a softer structure.³⁰ The hybrid PBE0 functional was shown to improve structural properties beyond that of PBE,³¹ and the addition of van der Waals to PBE0 provides an even greater increase in accuracy.³² Recently, the use of a dielectric-dependent hybrid (DDH) density functional was also shown to yield good agreement with experimental structural and electronic properties.³³ There is evidence that van der Waals functionals can be used to increase the accuracy of more complicated functionals such as the meta-GGA Strongly Constrained and Appropriately Normed (SCAN)³⁴ functional (SCAN+rVV10).³⁵ However, this is not necessarily improving all properties. Wiktor et al. found that it could decrease the accuracy of the structural properties of liquid water.³⁶ In the meta-GGA category, the SCAN functional³⁴ was recently shown to provide an accurate description of liquid water when used at a temperature of 330 K.^30,37 The choice of a simulation temperature of 330 K was made to approximately account for nuclear quantum effects (NQEs).^38,39 The importance of nuclear quantum effects (NQEs) in the description of water was discussed as early as 1985.^40,41 The ad hoc approach of increasing simulation temperature to account for NQEs has been discussed by Ruiz Pestana et al.⁴² We adopt this approach in the present work despite its shortcomings, due to the high cost of more accurate treatments of NQEs, and in order to allow for direct comparisons with previous work. Various combinations of density functionals have also been proposed, notably mixing GGA, meta-GGA, or hybrid functionals^22,43 with van der Waals functionals. A recent liquid water simulation using PBE0 and van der Waals terms led to structural properties in good agreement with experiment.³² A large number of studies of liquid water using several different exchange-correlation functionals have been completed over time.^{19,32,37,44–46} We refer to the recent review by Gillan et al.¹⁸ for a detailed account. Despite the large body of results published on first-principles simulations of water, a quantitative measure of the uncertainty of the results due to the statistical nature of molecular dynamics simulations has rarely been provided.

In this paper, we focus on the evaluation of the accuracy of the SCAN functional for the description of the structural and electronic properties of liquid water. To this end, we adopt a similar systematic methodology recently employed to benchmark the accuracy of the PBE functional at 400 K (PBE400).¹⁹ This approach is based on performing multiple independent MD simulations in order to obtain estimates of the variability of computed physical quantities. We performed a total of 16 independent simulations, each having a duration of 43.5 ps (cumulative duration of 696 ps) using the SCAN functional at 330 K (the resulting dataset is hereafter referred to as SCAN330). We compare the results of these simulations to those of the PBE400 dataset and to experimental data, focusing on structural parameters, electronic polarizability, as well as infrared and Raman spectra. Although the PBE functional used at 400 K is not representative of the state of the art in DFT simulations of water, this choice is motivated by the fact that PBE has been used with an elevated temperature by several authors in the past in an attempt to reproduce experimental structural data. Furthermore, the availability of the large PBE400 dataset allows for a comparison that includes a measure of the statistical uncertainty of simulations. Our results show that the SCAN density functional used at a simulation temperature of 330 K provides an improved overall description of liquid water—for both structural and electronic properties—compared to the PBE functional used at 400 K.

A series of 16 independent first-principles molecular dynamics simulations were performed. Each simulation included 64 D₂O molecules, placed in a cubic unit cell of size a = 23.46 (a.u.), consistent with the experimentally known density of heavy water 1.11 g/cm³. The electronic structure calculations were performed using the Strongly Constrained and Appropriately Normed (SCAN)³⁴ exchange-correlation functional implemented in the Qbox first-principles molecular dynamics code.^47,48 We use a plane wave energy cutoff of 85 Ry for MD simulations and an increased cutoff of 120 Ry when computing the stress tensor from the selected snapshots. We use norm-conserving pseudopotentials generated using the procedure of Hamann, Schlüter, Chiang⁴⁹ and Vanderbilt⁵⁰ (HSCV). This set of pseudopotentials was generated using the PBE functional. We assume that the use of a PBE functional (rather than the SCAN functional) when generating the pseudopotentials causes a negligible error in our calculations. This assumption is consistent with the work of others using the SCAN functional.^30,35,37,42

The molecular dynamics simulations were performed within the Born-Oppenheimer approximation, with the electronic ground state recomputed at every MD step. Heavy water was used to allow for large MD time steps. All simulations were performed with a time step of 10 (a.u.) (≃0.24 fs), and 1 800 000 time steps were performed resulting in ≃43.5 ps per trajectory. The Bussi-Donadio-Parrinello (BDP) thermostat⁵¹ was used with a target temperature of 330 K and a relaxation time of 5000 (a.u.) (≃0.12 ps).

In order to minimize the amount of time needed to achieve thermal equilibrium, the simulations were started from the end configurations of 16 independent simulations taken from the PBE400 dataset.^19,52 Using previous PBE simulations as a starting point relies on the assumption that PBE configurations are close to SCAN configurations, therefore reducing the amount of simulation time that has to be discarded before analyzing the data. Our choice of 16 simulations of 43.5 ps each was motivated by the need to have sufficiently long simulation time to enable the calculation of vibrational properties. Using more simulations of shorter duration would also lead to discarding more simulation time, since the initial equilibration period would have to be discarded for each simulation. We analyzed the convergence of our simulations using two different methods. First, we examined the Kohn-Sham energy vs simulation time, as shown in Fig. 1. This figure suggests that after ≃2–3 ps of simulation, the system achieved a state representative of its long term behavior. In order to get a more quantitative estimate of the point of equilibration, we used the marginal standard error rule (MSER) heuristic as described by Preston White,⁵³ 

where Y [i] is the ith element of a dataset Y of size N, d is the number of data elements discarded from the dataset (starting from the beginning of the dataset), and Ȳ_d is the average of the data Y after discarding d elements. The MSER quantifies the change in variance as data from the beginning of the sequence is discarded. If a simulation begins with outlier data and eventually converges to a set range, then the MSER will begin with a high value, decrease, and then increase again. The minimum of the MSER value indicates the optimal point at which to discard data so as to minimize the variance of the dataset. As shown in Fig. 2, the MSER reaches a minimum after about 1 ps. Therefore, in combination with the previous (visual) analysis of the Kohm-Sham energy time evolution, we conclude that discarding about 2.5 ps of data would be safe. However, out of extra caution, we decided to discard the first 5 ps of data. In the rest of this paper, any computed property does not include the first 5 ps of data. As an additional measure of the convergence of our simulations, we have computed the autocorrelation functions of the Kohn-Sham energy (averaged over each 0.5 ps run) in all 16 samples. The autocorrelation functions decay within 2–3 ps to small values that are within the limits of statistical significance (see the supplementary material).

Our analysis of the MD simulations starts with the evaluation of the average number of hydrogen bonds per water molecule. We follow the definition of a hydrogen bond proposed by Luzar and Chandler.⁵⁴ We computed the number of hydrogen bonds, averaged over the 64 molecules, during the last 38.7 ps (1 600 000 steps) of simulation for all trajectories. We found that the SCAN330 simulations had an average of 3.621 with a standard deviation of 0.11, compared to PBE400 which had an average of 3.566 and a standard deviation of 0.12.¹⁹ This indicates that using the SCAN functional at 330 K leads to an overall slightly stronger intermolecular bonding strength than the PBE functional at 400 K. This can be attributed both to the difference in simulation temperature and to the tendency of the SCAN functional to form stronger hydrogen bonds.

To further explore the structural properties of the simulated liquid, we computed atomic pair correlation functions. Figure 3 shows the oxygen-oxygen pair correlation function g_OO(r) for each of the 16 SCAN330 trajectories (blue) along with the corresponding values from the 32 PBE400 trajectories (red), and compared with experimental data (black). Figure 4 shows the same data averaged over all trajectories. The error estimate in the experimental data given by Soper⁵⁵ is shown as a green shaded region. The pair correlation functions were computed by binning the data into intervals of 0.05 (a.u.). When comparing computed data with experimental data, we should consider the fact that nuclear quantum effects (approximated in our calculations by a modified temperature) might further affect the computed values. With this caveat in mind, it can be seen that PBE400 results are closer to experimental data than the SCAN330 results at the first g_OO(r) maximum, whereas SCAN330 is more closely aligned with experimental data at the first minimum and at the second maximum. The availability of independent trajectories provides an estimate of the variability that should be expected in the g_OO(r) function when using simulations of approximately 38.7 ps and 64 molecules. The O–H pair correlation functions are shown in Fig. 5, which displays the expected large intramolecular O–H bond contribution. Figure 6 shows the same data averaged over all trajectories. While, in general, PBE400 is slightly closer to the experimental curve than SCAN330, neither SCAN330 nor PBE400 matches the second maximum very well in terms of magnitude, although they both match the second minimum and third maximum reasonably. Such disagreement with the experiment is expected for the g_OH(r) function since nuclear quantum effects are known to have a large impact on this function. A proper inclusion of NQEs using, e.g., path-integral Monte Carlo simulations rather than MD simulations is expected to reduce the amplitude of the second and third peaks in g_OH(r), as was demonstrated, e.g., by Medders et al. using the MB-pol force field.⁹ Finally, the H–H pair correlation is shown in Fig. 7 and averaged data over all trajectories in Fig. 8. Similar to the case of the O–H pair correlation, a large first peak describes the strong intramolecular correlation. Both SCAN330 and PBE400 yield similar results. As in the case of the g_OH(r) function, the magnitude of the second maximum is expected to be reduced by a proper treatment of the NQE.⁹ We note that both SCAN330 and PBE400 slightly underestimate the position of the second maximum. In addition to pair correlation functions, we also examined the O–O–O bond angle distribution shown in Fig. 9. The same data averaged over all trajectories is shown in Fig. 10. The angle formed by triplets of oxygen atoms was computed for atoms included in spheres centered on an oxygen atom with a radius of 3.2 Å (SCAN330) and 3.25 Å (PBE400). The sphere radii were chosen so as to have an average oxygen coordination number of 4. Comparisons with the experimental bond angle distribution are limited by the fact that experimental data are inferred using a model and not measured directly.⁵⁷ With this caveat in mind, it appears that SCAN330 yields an overall improvement over PBE400. A close examination shows that one SCAN330 trajectory represents an outlier and falls in the middle of the PBE400 data. This demonstrates the importance of using multiple independent trajectories, since incorrect conclusions might be drawn on the basis of a single 39 ps trajectory. Based on theanalysis of pair correlation functions, it appears that the SCAN functional provides a clear improvement over the PBE functional since it allows for simulations at a much more realistic temperature in order to approximate experimental data at ambient conditions. Overall, our computed structural properties are in agreement with the work of Chen et al.³⁷ and Zheng et al.,³⁰ but our results provide a measure of the variability of the data. For comparison, Ref. 37 includes a single simulation of 30 ps and Ref. 30 includes a single simulation of 50 ps, while the present study includes 16 separate simulations of 43.5 ps each (696 ps in total). Both the work of Chen and our own study discard the first 5 ps of data for equilibration, while Ref. 30 discards 8 ps. Beyond the analysis of pair correlation functions, a more detailed study of the structure of liquid water has been proposed by Gasparotto et al.⁵⁸ who discussed topological coordination defects and their correlations in ab initio simulations at various temperatures. A similar study on SCAN simulations would be of interest and will be considered in future work.

The ability of a density functional to reproduce the correct equilibrium density of liquid water at ambient conditions is a very useful criterion to identify accurate density functionals. This feature is critical in order to enable simulations of complex systems involving water, e.g., solvated ions or wet surfaces, for which the determination of the simulation volume cannot simply be inferred from experimental data but must rely on NpT simulations. The well-documented^18,31,59 property of the PBE functional (used at T = 400 K) to underestimate the density by about 15% has been a major obstacle to its use for solvated systems. Adding van der Waals corrected density functionals has been shown to yield equilibrium densities close to the experimental value.⁶⁰ The evaluation of the equilibrium density can, in principle, be carried out by performing an NpT simulation at ambient pressure and by measuring the resulting average volume. This procedure was followed by Chen et al.³⁷ who reported an average density of 1.05 g/cm³ for the SCAN functional. However, NpT simulations including only 64 molecules are subject to large volume fluctuations that result in large uncertainties in the equilibrium volume. We chose instead to compute the average stress tensor in a unit cell of fixed size (NVT) using the experimental volume. A comparison of the average stress allows one to discuss the relative merits of density functionals, using as a reference the fact that the PBE functional at 400 K yields an overestimate of the volume by about 15%.^37,46,59,61 We performed an analysis of the stress tensor on each trajectory by sampling every 2.5 ps starting from 5 ps for the SCAN330 dataset and 15 ps for the PBE400 dataset. The stress tensor was computed using an increased plane wave energy cutoff of 120 Ry. A histogram of the computed average of the diagonal elements of the stress tensor is shown in Fig. 11. The PBE and SCAN functionals both yield bell-shaped distributions centered near zero. Note that an ideal exchange-correlation functional should yield a distribution centered at zero. With an average of −0.29 GPa and a standard deviation of 0.43 GPa, SCAN330 is closer to zero than PBE400 which has an average stress of 0.54 GPa with a standard deviation of 0.44 GPa. This indicates that SCAN330 provides an improved description of the liquid density, although with an error of the opposite sign to PBE400, and a residual stress of about half that of PBE400. Assuming a similar bulk modulus of 2.2 GPa for SCAN330 and PBE400, we infer an equilibrium density of approximately 1.08 g/cm³ for the SCAN functional at 330 K.

The computation of the diffusion coefficient of water is known to be extremely sensitive to simulation parameters and, in particular, to finite-size effects.⁶² The high cost of DFT simulations requires the use of small unit cells, leading to a large uncertainty in the computed diffusion coefficient.¹⁹ We have computed the mean square displacement (MSD) of oxygen atoms vs time for the 16 SCAN330 trajectories. Figure 12 shows the MSD curve for each trajectory, along with the average over all trajectories. Unlike structural properties, MSD curves show a large variation when comparing one trajectory to another. As such, any comparison of the accuracy of the SCAN and PBE diffusive properties should be done cautiously. We computed the diffusion coefficient D_avg for all trajectories using a linear fit to the data after discarding the first 5 ps. For SCAN330, we find the average value of D_avg to be 1.26 × 10⁻⁵ cm²/s with a standard deviation of σ_D = 0.4 × 10⁻⁵. The standard deviation of the mean is $σDavg$ = σ_D/$N$⁠, where N is the number of samples, yielding $σDavg$ = 0.1 × 10⁻⁵. For PBE400, the values were computed by Dawson and Gygi¹⁹ who reported D_avg = 1.96 × 10⁻⁵ cm²/s, σ_D = 0.5 × 10⁻⁵, and $σDavg$ = 0.1 × 10⁻⁵cm²/s.

Finite-size effects on the diffusion coefficient of water were investigated by Yeh and Hummer⁶³ who used a TIP3P model to compute the diffusion coefficient in samples of sizes ranging from 128 to 2048 molecules. These authors used an additive correction proposed by Dünweg and Kremer⁶⁴ that depends on the computed viscosity of the system and allows for an extrapolation to samples of infinite size. An accurate computation of the viscosity of water would require exceedingly long first-principles simulations. We note that the computed viscosity used in Ref. 64 may differ substantially from the experimental value. Thus, the use of the experimental viscosity in the correction may introduce a large uncertainty. Furthermore, a comparison of the diffusion coefficient with the experimental value is made difficult by the strong temperature dependence of this quantity. The experimental value⁶⁵ varies from 1.15 × 10⁻⁵ cm²/s at 274 K to 3.58 × 10⁻⁵ cm²/s at 318 K. Considering the ad hoc correction to the effective simulation temperature of 30 K adopted in this work to account for nuclear quantum effects, it is unclear how a computed value of D—even corrected for finite-size effects—should be compared with the experiment. We therefore refrain from comparing our result directly with experiment, although we note that finite-size corrections would generally increase the value of the diffusion coefficient. Comparisons between different density functionals remain, however, valid using uncorrected values, assuming that simulations include the same number of molecules. The comparison of the SCAN330 and the PBE400 results indicate that SCAN330 yields slower diffusion than PBE400.

Beyond structural and dynamical properties, it is highly desirable to obtain an accurate description of the electronic properties of water, since this is likely to affect several properties including solvation of ions and wetting of surfaces. We have tested the ability of the SCAN and PBE functional to accurately describe electronic properties by computing the high-frequency dielectric constant (electric permittivity) as well as the infrared (IR) and Raman spectra. We computed the static electronic polarizability tensor α_ij on simulation snapshots separated by 2.5 ps. The tensor is computed by solving the Kohn-Sham equations in the presence of a constant total electric field applied to the sample in each of the x, y, and z directions and using a centered 3-point finite difference formula to compute the elements of α_ij according to

where d_i is the i component of the total dipole moment and E_j is the j component of the total electric field. The constant electric field was implemented following the electric enthalpy formulation of Souza and Vanderbilt.⁶⁶ The total dipole moment in the unit cell was computed using Maximally Localized Wannier functions (MLWF)^67,68 with the refinement scheme of Stengel and Spaldin.⁶⁹ The dielectric constant was then calculated as

where Ω is the unit cell volume and

The trajectories were sampled every 2.5 ps, after discarding the first 5 ps of simulation for SCAN330 (total of 272 samples) and discarding 15 ps for PBE400 (total of 576 samples). Figure 13 shows histograms of the high-frequency dielectric constant, compared to the experimental value of 1.77 reported in Ref. 70. It clearly shows that SCAN330 is significantly more accurate than PBE400 for the description of ϵ_∞. A comparison with the experimental value shows that the error is reduced by half when switching from PBE400 to SCAN330. We also computed the molecular dipole moments from the positions of the Wannier centers and compared them in the PBE400 and SCAN330 calculations. The average molecular dipoles are 3.04 D and 2.96 D for PBE and SCAN, respectively. This result is consistent with the findings of Chen et al.³⁷ who found a value of 2.97 D for the SCAN functional. A histogram of the molecular dipole values in PBE and SCAN shows that the two distributions have a similar width, but are shifted by approximately 0.1 D (see the supplementary material). The positions of the Wannier centers are also slightly different in the PBE400 and in the SCAN330 datasets. We have computed the distribution of the distance d_OW between the oxygen atom of each molecule and the nearest four Wannier centers. A comparison of this distribution in the PBE400 and SCAN330 datasets shows the characteristic double peak of the distribution corresponding to the centers located in the lone pair directions (d_OW ≃ 0.32 Å) and the centers located along the OH bond directions (d_OW ≃ 0.5 Å). The peaks of the PBE400 distribution are slightly shifted to larger distances by approximately 0.01 Å. Similarly, the spread of the Wannier functions is slightly shifted to larger values in the PBE400 dataset (see the supplementary material).

We completed our analysis of electronic properties by computing the IR and Raman spectra. The IR spectrum was computed from the Fourier transform of the electric dipole moment autocorrelation function sampled every 2.5 fs. We used the following definition of the IR absorption coefficient A(ω):

where n(ω) is the refractive index, V is the volume, β = 1/k_BT, and $M→$ is the total dipole moment. The computed infrared absorption spectrum obtained using Eq. (5) was filtered with a Gaussian kernel having a half width at half maximum (HWHM) of 36 cm⁻¹. The Raman spectrum was computed from the isotropic electric polarizability (α_iso) evaluated every 2.5 fs. The isotropic polarizability was obtained from the polarizability tensor computed by finite differences as described above. Due to the high cost of the computation of the polarizability, we used only the last 20 ps of each trajectory to compute the IR and Raman spectra and limited ourselves to only 8 trajectories for each functional. Figure 14 shows the IR absorption spectra for SCAN330, PBE400, and experimental data from Max and Chapados,⁷¹ and Fig. 15 shows the same spectra averaged over all trajectories. The experimental spectrum is arbitrarily scaled, while the computed spectra are shown on the same scale (see the supplementary material). The infrared spectrum of water has been previously computed by Xu et al.⁷² using the SCAN functional in a single simulation of 50 ps, using Car-Parrinello dynamics. Our results agree with those of Xu et al. and add a measure of the uncertainty affecting comparisons of the infrared spectrum. Focusing on the O-D stretch peak, we observe that SCAN appears to predict the location of the peak with good accuracy. The lower frequency modes are well modeled by both SCAN and PBE with neither being clearly better than the other. We note, however, that the agreement with the experimental peak positions would be further modified if nuclear quantum effects were properly taken into account. Previous studies^73–76 have shown that the inclusion of NQEs causes a substantial red shift of the O-D stretch mode. The observed blue-shift of this mode upon moving from PBE to SCAN therefore represents a partial improvement, which would, however, be red shifted again by the inclusion of NQEs. The improvement in peak positions for the SCAN functional can be traced to an improved description of the O-D stretch frequencies in the water monomer. A simple analysis of a single D₂O molecule was done in order to examine and compare the vibrational frequencies of water for SCAN and PBE, and to validate our choice of plane wave cutoff. The summary is shown in Table I. There is no meaningful difference in the three PBE frequencies for cutoffs ranging from 85 to 200 Ry. Tests confirmed that SCAN and PBE have the same convergence trend, and therefore, 85 Ry is sufficient for our MD simulations. We find that SCAN provides better accuracy than PBE for all three vibrational modes. Interestingly, PBE underestimates all three frequencies, while SCAN underestimates two and overestimates the third.

The failure of the SCAN functional to correctly reproduce the high-frequency part of the spectrum (when considering NQEs) could be attributed to the overestimation of the hydrogen bond strength observed, e.g., in calculations of the binding energy of weakly interacting closed shell systems.⁷⁸ 

The Raman intensity was computed by evaluating the system’s polarizability at intervals of 2.5 fs for eight independent trajectories for each choice of density functional. Polarizabilities were computed using a finite constant electric field in the x, y, and z directions and using a finite-difference expression of the polarizability. The Raman intensity was computed following the approach of Wan et al.⁷⁹ Figure 16 shows the Raman spectra for SCAN330, PBE400, and experimental data from Pattenaude et al.,⁸⁰ and Fig. 17 shows the same spectra averaged over all trajectories. The experimental data are arbitrarily scaled, while the computed spectra are shown on the same scale (see the supplementary material). The SCAN functional predicts the location of the O-D stretch peak accurately, although the inclusion of nuclear quantum effects will likely affect the peak position as in the case of the IR spectrum. Our results for the PBE400 dataset are in good agreement with the results reported by Wan et al.⁷⁹ As in the case of the IR spectrum, we find that the SCAN functional gives an improved description of the Raman spectrum when compared to PBE results.

Throughout our analysis, we found that the SCAN meta-GGA functional provides a good description of the structural and electronic properties of liquid water at a moderate computational cost. When compared to the simple approach of using PBE at 400K, SCAN brings significant improvements in all computed properties. While it is difficult to clearly separate the parts of electronic correlations that are included in either functional, SCAN as a meta-GGA has additional flexibility to satisfy constraints in reproducing known exact results. In particular, the inclusion of some of the intermediate-range van der Waals interaction in SCAN is credited by its authors to yield substantially improved binding energies in the S22 dataset of weakly interacting closed-shell molecules.⁸¹ The observation that an increase in the simulation temperature in PBE is needed to compensate for the lack of vdW interactions in liquid water is consistent with similar findings by Gasparotto et al.⁵⁸ who discussed the same effect for the BLYP functional. The semilocal character of SCAN, however, makes it unable to properly describe long-range van der Waals interactions, which are necessary for a more accurate description of liquid water. Such long-range interactions are included in high-level classical models such as the MB-pol force field, which is derived from high-accuracy quantum chemical calculations.⁹ The calculated infrared and Raman spectra show that SCAN, although it brings an improvement over PBE, does not reach a good agreement with the experiment if one accounts for the fact that nuclear quantum effects would further shift the high-frequency part of the spectra to lower frequencies. The SCAN functional is also found to yield an improved description of the equilibrium volume when compared to PBE. Our results also confirm the need for a better treatment of the nuclear quantum effects in water beyond the ad hoc approach of increasing the simulation temperature. Overall, the substantial variability of results obtained with 64 molecules and moderate simulation times (e.g., 40 ps) demonstrate that ensemble simulations are important in order to make meaningful comparisons between density functionals. We provide all simulation data and trajectories at http://quantum-simulation.org.

See the supplementary material for the infrared and Raman spectra normalized by the amplitude of the largest peak, the autocorrelation functions of the Kohn-Sham energy, the distribution of molecular dipole moments, the distribution of the oxygen-Wannier center distances, and the distribution of the Wannier function spreads.

We thank Jianwei Sun for providing the original implementation of the SCAN exchange-correlation functional for verification of our implementation, the authors of Ref. 71 for providing a digital version of the experimental infrared absorption spectrum, and G. Galli for a critical reading of the manuscript. This work was supported by the Midwest Integrated Center for Computational Materials (MICCoM) as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division (Grant No. 5J-30161-0010A). This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357. One of us (F.G.) gratefully acknowledges the hospitality of the Pritzker School of Molecular Engineering at the University of Chicago, as well as the Laboratoire de Théorie et Simulation des Matériaux (THEOS) and NCCR Marvel at the Ecole Polytechnique Fédérale de Lausanne, Switzerland, where part of this work was completed.