# Ab initio thermodynamics of liquid and solid water

^{a}Laboratory of Computational Science and Modeling, Institute of Materials, École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland;^{b}Universität Göttingen, Institut für Physikalische Chemie, Theoretische Chemie, 37077 Göttingen, Germany;^{c}International Center for Advanced Studies of Energy Conversion, Universität Göttingen, 37073 Göttingen, Germany;^{d}Faculty of Physics, University of Vienna, 1090 Vienna, Austria

See allHide authors and affiliations

Edited by Pablo G. Debenedetti, Princeton University, Princeton, NJ, and approved December 3, 2018 (received for review September 4, 2018)

## Significance

A central goal of computational physics and chemistry is to predict material properties by using first-principles methods based on the fundamental laws of quantum mechanics. However, the high computational costs of these methods typically prevent rigorous predictions of macroscopic quantities at finite temperatures, such as heat capacity, density, and chemical potential. Here, we enable such predictions by marrying advanced free-energy methods with data-driven machine-learning interatomic potentials. We show that, for the ubiquitous and technologically essential system of water, a first-principles thermodynamic description not only leads to excellent agreement with experiments, but also reveals the crucial role of nuclear quantum fluctuations in modulating the thermodynamic stabilities of different phases of water.

## Abstract

Thermodynamic properties of liquid water as well as hexagonal (Ih) and cubic (Ic) ice are predicted based on density functional theory at the hybrid-functional level, rigorously taking into account quantum nuclear motion, anharmonic fluctuations, and proton disorder. This is made possible by combining advanced free-energy methods and state-of-the-art machine-learning techniques. The ab initio description leads to structural properties in excellent agreement with experiments and reliable estimates of the melting points of light and heavy water. We observe that nuclear-quantum effects contribute a crucial

- ab initio thermodynamics
- machine-learning potential
- water
- density functional theory
- nuclear quantum effects

Liquid water and ice are ubiquitous on Earth, and their thermodynamic properties have important consequences in the climate system (1), the ocean, biological cells (2), refrigeration, and transportation systems. The solid phase that is stable at ambient pressure is ice Ih, whose hexagonal crystal structure is reflected in the sixfold symmetry of snowflakes. The cubic form, Ic, is a metastable ice phase whose relative stability with respect to ice Ih plays a central role in ice cloud formation in the Earth’s atmosphere (3⇓–5) but is extremely difficult to measure experimentally (1).

Despite the simple chemical formula, *i*) the shortcomings of common water models including conventional force fields (6) and (semi-)local density functional theory (DFT) approaches (7⇓–9), (*ii*) proton disorder in ice, and (*iii*) the importance of nuclear quantum effects (NQEs) (10). In particular, calculating the chemical potential difference

Water and ice have been described with varying success by invoking approximations of differing severity, including simple electrostatic dipole models for the energetics of proton ordering (14), force-field-based path-integral molecular dynamics (PIMD) studies (15⇓⇓–18), first-principles quasiharmonic (QHA) (17, 19), and vibrational self-consistent field (VSCF) (13, 20) studies which provide an approximate upper bound for

In this study, we make theoretical predictions of thermodynamic properties of ice and liquid water at a hybrid DFT level of theory, taking into account NQEs, proton disorder, and anharmonicity. This is made possible by exploiting advances in machine-learning (ML) techniques to avoid the prohibitively large computational expenses otherwise incurred by extensively sampling phase space by using first-principles methods. In particular, we use sophisticated thermodynamic integration (TI) techniques to accurately and rigorously compute the chemical potential difference between ice Ic and Ih and between ice Ih and liquid water.

## First-Principles Thermodynamics

As the underlying electronic structure description, we use the hybrid revPBE0 (21⇓–23) functional with a Grimme D3 dispersion correction (24, 25), which has been demonstrated to accurately predict the structure, dynamics, and spectroscopy of liquid water in molecular dynamics (MD) and PIMD simulations (26). revPBE0-D3 predicts a difference in lattice energy between the most stable proton-ordered forms of ice Ic and Ih of *SI Appendix* for further details), which is consistent with diffusion Monte Carlo predictions of

Since thorough sampling of the phase space of water at the revPBE0-D3 level of theory is prohibitively expensive, we sample the phase space using a surrogate ML PES, **1** is rendered particularly affordable and robust by the high fidelity of our surrogate ML PES, which substantially exceeds that obtained from empirical force fields or local DFT calculations, which were previously used as implicit surrogates (28, 29). Eq. **1** is the central formula of our approach: It not only enables accurate and efficient free-energy estimation at the ab initio level by delegating phase-space sampling to cheap surrogate models, but also provides a general way for benchmarking and calibrating the ML potentials.

## Neural Network Potential for Water

We constructed a flexible and fully dissociable neural network (NN) potential for bulk liquid water and ice following the framework of Behler and Parrinello (30⇓–32) using the RuNNer code (33), which was trained on the basis of revPBE0-D3 energies and forces for 1,593 diverse reference structures of 64 molecules of liquid water computed by using the CP2K package (34). Further details regarding the DFT calculations, comparison with the energies computed by using VASP (35), and the training and validation of the NN potential can be found in *SI Appendix*. We have released this NN potential along with its training set (*SI Appendix*, Datasets S1 and S2).

The revPBE0-D3-based NN potential describes the density (Fig. 1) and structural properties of water (Fig. 2) in very good agreement with experiments. Fig. 1 shows density isobars computed for ice Ic, ice Ih, and liquid water considering both the case of classical and quantum-mechanical nuclei. Fig. 1 highlights that (*i*) the predicted densities of liquid water and ice Ih and Ic agree with experiment to within 3%; (*ii*) the predicted thermal expansion coefficients show excellent agreement with experimental data; and (*iii*) the temperature of maximum density for liquid water matches the experimental value of

Fig. 2, *Top*, shows that NQEs have a slight destructuring effect on the oxygen–oxygen (O–O) radial distribution function (RDF), bringing it in excellent agreement with experiments from X-ray diffraction measurements (37), as also seen in the ab initio (PI)MD calculations with revPBE0-D3 (26). This destructuring has been observed in simulations using other DFT functionals (42) as well as empirical water models (43, 44) and was rationalized as a result of competing quantum effects (16, 45). Fig. 2, *Middle* and *Bottom*, further shows that NQEs cause a significant broadening of the oxygen–hydrogen (O–H) and hydrogen–hydrogen (H–H) RDFs, especially around their respective first peaks, which plays a predominant role in ensuring the match between the simulations and experiment. It is worth noting that the agreement between the NN and the experimental RDFs in Fig. 2 is significantly better compared with most benchmarked empirical water models and DFT functionals (46, 47).

## Promoting ML Potential to DFT

Despite the excellent performance of the NN potential, the fitting strategy, the finite cutoff radii applied to the description of atomic environments, and possible “holes” in the training set (48) inevitably lead to small residual errors compared with the underlying first-principles reference. To assess their significance, we have trained a collection of NN potentials using different training sets and/or initial random seeds, which demonstrates that predictions of the chemical potential difference between ice Ic and Ih from two different NN potentials can be as large as *SI Appendix*, Fig. S4 for further detail). Promoting the results to the DFT level eliminates these residual errors and any dependence on the specific NN potential used. This allows us to achieve submillielectronvolt accuracy in free energies (as required to resolve the greater stability of ice Ih compared with Ic) and to make unbiased properties predictions at the reference ab initio level of theory in general.

The temperature-dependent DFT corrections to the NN chemical potentials of different phases of water, **1**) performed on 64-molecule systems, are shown in Fig. 3. For each ice phase (Ic and Ih), 16 different proton-disordered initial configurations with zero net polarization are generated by using the Hydrogen-Disordered Ice Generator (49). The SD of the potential energy for the 16 proton-disordered ice Ic configurations is *i*) the proton order is effectively “frozen-in” at the timescales available to simulation (50) and (*ii*) there are significant differences between

## Results and Discussion

### The Relative Stability of Hexagonal and Cubic Ice.

We follow the workflow in *Materials and Methods* (see Fig. 7 for illustration) to evaluate the chemical potential difference

Note that the speed and linear scaling of the NN potential allow us to simulate ice systems containing as many as 768 water molecules. Such system size is not only essential to represent the wide spectrum of possible local arrangements realized in proton-disordered ice, but also important for averaging over different proton-disordered structures when correcting for the chemical differences between the NN potential and revPBE0-D3, as demonstrated by the spread of

NQEs are taken into account by integrating the quantum centroid virial kinetic energy *Materials and Methods* and see Fig. 6). We perform NN-based PIMD simulations within the NPT ensemble and assess the impact of NQEs on the chemical potential at the NN level using

Fig. 4 shows that the NN predictions of *SI Appendix* for further details). Assuming that the heat of transition from ice Ic to ice Ih is approximately constant over the temperature range 200–300 K, the temperature dependence of **2**) a transition enthalpy of

### The Relative Stability of Hexagonal Ice and Liquid Water.

Now, we compute the difference in chemical potential

We first compute *SI Appendix*, Fig. S3). Afterward, the calibration terms for chemical potentials

Fig. 5 shows

NQEs lower the melting point of

## Conclusions

We show that a revPBE0-D3 description of the electronic structure predicts properties for ice Ih, ice Ic, and liquid water that are in excellent quantitative agreement with experiment. This is made possible by using a ML potential as an intermediate surrogate model and by using advanced free-energy techniques. We not only rigorously compute but also quantitatively analyze the individual contributions from NQEs, proton disorder, and anharmonicity.

This study demonstrates that it is possible to achieve a submillielectronvolt level of statistical accuracy in predicting the thermodynamic properties of a complex system such as water at a hybrid DFT level of theory. The idea of using ML potentials as sampling devices significantly broadens the applicability and prowess of electronic-structure approaches, making it affordable to use them in the accurate computations of free energies and other thermodynamic properties. The overall framework and the free-energy methods described here provide a general, accurate, and robust way for first-principles predictions of thermodynamic properties of a plethora of physical systems, such as pharmaceutical compounds, hydrogen storage materials, hydrocarbons, and metallic alloys.

## Materials and Methods

### Simulation Details.

The density isobar in Fig. 1 is computed by using both classical MD and PIMD simulations in the NPT ensemble for ice Ic, ice Ih, and liquid water systems of 64 molecules. We confirm that the equilibrium density computed with 64 water molecules in classical MD simulations is consistent with the values obtained for systems with ∼2,000 molecules. All MD simulations and PIMD simulations that use 56 beads are performed by using the i-PI code (60) in conjunction with LAMMPS (61) with a NN potential implementation (62, 63).

### Interface Pinning.

The interface pinning simulations (56) are performed by using the PLUMED code (64) on an ice–liquid system containing 5,760 molecules at temperatures ranging from 250 to 300 K and pressure 1 bar, using the NN potential.

### Accounting for NQEs.

NQEs on the chemical potential difference between ice Ic and ice Ih at 200, 250, 273, and 300 K are taken into account by integrating the quantum centroid virial kinetic energy

### Workflow for Computing Δ μ Ih → Ic .

Here, we describe the workflow for computing absolute Gibbs free energy and thereby the chemical potential of an ice system. The first step is a TI from a harmonic reference to a classical ice system (**1** is included. Finally, to describe ab initio ice with quantum-mechanical nuclei (**3**). As an alternative strategy, one can also follow the TI route **1**, which is more costly.

### Datasets.

The NN potential for water based on revPBE0-D3, the training set for the potential, and all necessary simulation input files are included in *SI Appendix* and are available at https://archive.materialscloud.org/2018.0020/v1.

## Acknowledgments

We thank Matti Hellström for providing the VASP data and Gerit Brandenburg for benchmarking the revPBE0+D3 functional. B.C. was supported by Swiss National Science Foundation Project 200021-159896; a partial funding of a visit in Göttingen by the International Center for Advanced Studies of Energy Conversion; and generous allocation of CPU time by the Swiss National Supercomputing Center under Project s787. M.C. and E.A.E. were supported by the European Research Council under European Union’s Horizon 2020 Research and Innovation Programme Grant 677013-HBMAP. J.B. was supported by DFG Heisenberg Professorship BE3264/11-2. C.D. was supported by Austrian Science Fund Spezialforschungsbereich Vienna Computational Materials Laboratory Project F41.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: bingqing.cheng{at}epfl.ch.

Author contributions: B.C., J.B., C.D., and M.C. designed research; B.C. performed research; B.C. and E.A.E. analyzed data; and B.C., E.A.E., J.B., C.D., and M.C. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: The NN potential has been deposited on GitHub and is available at https://github.com/BingqingCheng/ab-initio-thermodynamics-of-water.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1815117116/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- ↵
- ↵
- ↵
- Lupi L, et al.

- ↵
- Kuhs WF,
- Sippel C,
- Falenty A,
- Hansen TC

_{c}” Proc Natl Acad Sci USA 109:21259–21264. - ↵
- ↵
- Morales MA, et al.

- ↵
- Zhang C,
- Wu J,
- Galli G,
- Gygi F

- ↵
- ↵
- Markland TE,
- Ceriotti M

- ↵
- Herrero CP,
- Ramírez R

- ↵
- ↵
- Engel EA,
- Monserrat B,
- Needs RJ

- ↵
- Lekner J

- ↵
- ↵
- ↵
- ↵
- Cheng B,
- Behler J,
- Ceriotti M

- ↵
- ↵
- Engel EA,
- Li Y,
- Needs RJ

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Marsalek O,
- Markland TE

- ↵
- Macher M,
- Klimeš J,
- Franchini C,
- Kresse G

- ↵
- Grabowski B,
- Ismer L,
- Hickel T,
- Neugebauer J

- ↵
- Glensk A,
- Grabowski B,
- Hickel T,
- Neugebauer J

- ↵
- ↵
- Behler J

- ↵
- Morawietz T,
- Singraber A,
- Dellago C,
- Behler J

- ↵
- Behler J

- ↵
- Lippert G,
- Hutter J,
- Parrinello M

- ↵
- Kresse G

- ↵
- ↵
- Skinner LB,
- Benmore C,
- Neuefeind JC,
- Parise JB

- ↵
- ↵
- Chen W,
- Ambrosio F,
- Miceli G,
- Pasquarello A

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Cisneros GA, et al.

- ↵
- Behler J

- ↵
- Matsumoto M,
- Yagasaki T,
- Tanaka H

- ↵
- Drechsel-Grau C,
- Marx D

- ↵
- Cheng B,
- Ceriotti M

- ↵
- Molinero V,
- Moore EB

- ↵
- Cheng B,
- Dellago C,
- Ceriotti M

- ↵
- Reddy SK, et al.

- ↵
- Carr THG,
- Shephard JJ,
- Salzmann CG

- ↵
- Pedersen UR,
- Hummel F,
- Dellago C

- ↵
- Haji-Akbari A,
- Debenedetti PG

- ↵
- ↵
- Rossi M,
- Fang W,
- Michaelides A

- ↵
- Kapil V, et al.

*Computer Phys Commun*, 10.1016/j.cpc.2018.09.020. - ↵
- Plimpton S

- ↵
- Singraber A,
- Behler J,
- Dellago C

- ↵
- Singraber A

*CompPhysVienna/n2p2: Neural Network Potential Package 1.0.0*(University of Vienna, Vienna). - ↵
- ↵
- ↵
- Cheng B,
- Ceriotti M

- ↵
- Cheng B,
- Paxton AT,
- Ceriotti M

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences