Extracting ice phases from liquid water: why a machine-learning water model generalizes so well
Bartomeu Monserrat, Jan Gerit Brandenburg, Edgar A. Engel, Bingqing Cheng
EExtracting ice phases from liquid water:why a machine-learning water model generalizes so well
Bartomeu Monserrat,
1, 2
Jan Gerit Brandenburg,
3, 4
Edgar A. Engel, and Bingqing Cheng
5, 2, ∗ Department of Materials Science and Metallurgy, University of Cambridge,27 Charles Babbage Road, Cambridge CB3 0FS, United Kingdom Cavendish Laboratory, University of Cambridge,J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom Interdisciplinary Center for Scientific Computing, University of Heidelberg,Im Neuenheimer Feld 205A, 69120 Heidelberg, Germany Chief Digital Organization, Merck KGaA, Frankfurter Str. 250, 64293 Darmstadt, Germany Department of Chemistry, University of Cambridge,Lensfield Road, Cambridge, CB2 1EW, United Kingdom (Dated: June 25, 2020)We investigate the structural similarities between liquid water and 53 ices, including 20 knowncrystalline phases. We base such similarity comparison on the local environments that consist ofatoms within a certain cutoff radius of a central atom. We reveal that liquid water explores the localenvironments of the diverse ice phases, by directly comparing the environments in these phases usinggeneral atomic descriptors, and also by demonstrating that a machine-learning potential trained onliquid water alone can predict the densities, the lattice energies, and vibrational properties of theices. The finding that the local environments characterising the different ice phases are found inwater sheds light on water phase behaviors, and rationalizes the transferability of water modelsbetween different phases.
Keywords: water, phase behavior, machine learning, polymorphism, atomic descriptors
I. INTRODUCTION
The number of experimentally observed and theoret-ically predicted phases of water seems to be ever grow-ing [1]. The most ubiquitous phase on Earth, liquidwater, has many intriguing properties, including a den-sity maximum at ◦ C and ambient pressure, volume ex-pansion upon freezing, unusually high surface tension,melting and boiling point [2]. Liquid water exhibits nolong-range order and its local structure is difficult toquantify and yet intricately related to its unique prop-erties [2–5]. Beside the liquid, the various ice phasesin the complex phase diagram of water are made fromdistinct local atomic environments [1], which lead toa large spread in their densities, lattice energies, andother thermodynamic as well as kinetic properties [1, 6].Apart from the direct connection with physical proper-ties, the local structures in water are also related to thetransition paths between the phases [7, 8].One intriguing question thus is the structural rela-tionship between the different ice phases, and betweenice and liquid water. This is not an easy topic to inves-tigate, however, due to the structural complexity of theliquid [4, 9, 10] and the large number of ice phases [1].In this work, we exploit the state-of-the-art advancesin machine learning (ML) for chemistry and materials, ∗ Correspondence email address: [email protected] in order to compare the local environments in variousphases of water in a general and systematic manner.More specifically, we first curate a dataset consisting of53 representative phases of ice including all the knownphases (see Sec V B), whose densities range from 0.7 to . / mL . Then we demonstrate that the local atomicenvironments found in liquid water cover the ones ob-served in all these ice phases, using a universal and au-tomated framework for comparing the local similarities.As a consequence of this inclusion, a machine-learningpotential (MLP) [11] that is only trained on liquid wa-ter accurately reproduces the ice properties includinglattice energies, mass densities, and phonon density ofstates. II. RESULTSA. Curated dataset of diverse water phases
We first select representative atomistic configurationsof diverse crystalline and liquid phases. We start from57 ice crystal structures, which include all the experi-mentally known ices. These were screened from an ex-tensive set of 15,859 hypothetical ice structures using ageneralised convex hull construction (an algorithm foridentifying promising experimental candidates) [12, 13].After rigorous geometry optimizations at zero pressure(see Sec. V C), we eliminated three defected phases andthe very high pressure phase X. Section V B describes a r X i v : . [ phy s i c s . c o m p - ph ] J un the dataset of the remaining 53 ice phases in more de-tail. Note that some structures (with particular hydro-gen arrangements) represent both a proton-ordered anda proton-disordered form: for example, one ice struc-ture prototypes both ice Ih and XI. We consider therespective minimum potential energy configurations ofthe ice phases, because they provide reasonable and re-producible approximations to the physical properties ofice, and serve as starting points for computing thermo-dynamic properties.Compiling a set of representative structures for liq-uid water is less straightforward, since the liquid per-sists over a wide range of temperatures and pressures.We consider 1,000 diverse 64-molecule snapshots of liq-uid water, which have previously served in training arecent MLP [11]. They were originally prepared usinga three-step process. Bulk liquid systems of 64 watermolecules were first equilibrated at high temperaturesand densities between 0.7 and . / mL . The resulting(de-correlated) configurations were then quenched usinga steepest decent optimization. Finally, the 1,000 moststructurally diverse structures were extracted from allthe collected liquid configurations using a farthest pointsampling algorithm. There are two reasons why the1,000 configurations are representative of liquid water.First, they were constructed in order to cover a largepart of the configurational space of possible atomic en-vironments in liquid water. Second, the MLP trainedusing these structures reproduces many properties ofwater very well, including the density isobar and radialdistribution functions at ambient pressure [11], whichmeans that the training set contains the necessary infor-mation for describing liquid water at ambient pressurein a data-driven manner. B. Direct comparison of the local environments
We employ the Smooth Overlap of Atomic Positions(SOAP) [14] local descriptors to represent the atomicenvironments (i.e. the displacements of all the neigh-bors within a cutoff radius r c around the central atom).Section V A provides more details regarding the rep-resentations. For each structure, we then compute itsglobal descriptors by taking the average of the local onesof all the atomic environments in that structure [15].As the global descriptors are high-dimensional, we usethe principal component analysis (PCA) to build a two-dimensional embedding to visualise the relative differ-ence (i.e. distances) between the structures. Essen-tially, a PCA map is linear projection that best pre-serves the variances of the high-dimensional Cartesiandistances of the dataset. Because only linear operationsare involved throughout, the local and the global de-scriptors can be meaningfully projected onto the samePCA map. We use these methodologies to analyse the 53 icephases and the 1,000 snapshots of the liquid. Fig. 1(a) shows the PCA map of the global descriptors ofall the structures: similar structures stay close on thismap, while distinct ones are farther apart. The hori-zontal principal axis is strongly correlated with density,suggesting that density variance is a dominant featureof the dataset. The ice phases and the liquid structuresare separated on the map, while ice structures that arecommonly considered to be similar (e.g. ice Ic and Ih)stay close together. The distinction between ices andliquids in the PCA map is to be expected consideringthe absence of long-range order in the latter.Fig. 1 (b) shows the projection of local descriptors ofall the atomic environments ( r c = 6 Å), onto the samePCA map. The local environments in liquids and icesare similar on this map. Furthermore, the crystallinenature of the ices leads to comparatively few distinctatomic environments, and they are almost completelycovered by the continuum realised in liquid water. Inother words, liquid water prototypes all atomic environ-ments pertinent to the 53 ice phases.
C. Predictions of the MLP on ices
The PCA maps provide a simple and general way ofcomparing and understanding the structural similari-ties, but choice of the representations and the lineardimensionality reduction inevitably lead to informationloss and distortion. As an alternative similarity com-parison, we explore how well a MLP [11] that is onlytrained on reference calculations for liquid water con-figurations describes diverse crystalline phases.This MLP is based on revPBE0 [17, 18] hybrid func-tional DFT calculations with the semi-classical D3 dis-persion correction [19]. The training set contains 1,593configurations: the first 1,000 are classical configura-tions as described above, the remaining 593 originatefrom path-integral molecular dynamics (PIMD) simula-tions at ambient conditions. We omit those PIMD con-figurations in the PCA analysis above because nuclearquantum effects [20] complicate the direct comparisonswith classical water, and also because those configura-tions had a very minor effect on the training of the MLPin our previous work [11]. The MLP uses an artificialneural network constructed according to the frameworkof Behler and Parrinello [21]. The total energy of thesystem is expressed as the sum of the individual contri-butions from the atom-centered environments of radius6 Å.Crucially, the success of the MLP hinges on the no-tion of “nearsightedness”: energy and forces associatedwith a central atom are largely determined by its neigh-bors, and the long-range interactions can be approxi-mated in a mean-field manner without explicitly consid-
Figure 1. PCA maps for the 53 ice phases and the 1,000 liquid water configurations. The geometries of the ice phaseshave been optimized by HSE-3c (see Sec. V C in Methods). If known, each ice structure is labelled by the name of theproton ordered and disordered phase, and otherwise by a number. Upper panel: each dot indicates a structure, which isdescribed by global descriptors. Lower panel: each small dot indicates the projection of a oxygen or hydrogen-centeredlocal environment. ering the far-away atoms. This notion underlies manyatomic and molecular force-field as well as most com-mon MLPs [22]. From this point of view, to capturethe energetics and dynamics of a phase of water, thekey is to predict the local atom-centered contributionsto the total energy the forces. In practice, this meansthe training set of the MLP needs to contain the essen-tial local atomic environments of the particular phase.Following this logic, we postulate that, if the liquid wa-ter contains all the local environments of the ice phases, a MLP trained exclusively on snapshots of liquid watershould also be able to describe the ice phases.To verify our hypothesis, we benchmark the perfor-mance of the MLP against reference DFT calculationsand experimental results (see more details in Meth-ods). The DFT references comprise (i) revPBE0-D3using cp2k with similar numerical settings as the cal-culations performed to generate the training referenceof the MLP, (ii) revPBE0-D3 using vasp and convergednumerical settings, and (iii) the hybrid HSE-3c as im-
Figure 2. A comparison between the lattice and energy of the 53 ice phases computed using the MLP based on revPBE0-D3, revPBE0-D3 DFT calculations employing vasp , revPBE0-D3 DFT calculations employing cp2k , and HSE-3c DFTcalculations using
CRYSTAL17 . The experimental references are taken from Ref. [16].Figure 3. A comparison between the phonon density ofstates (DOS) for the ice Ih (XI) phase computed usingrevPBE0-D3 DFT calculations employing vasp for a 8molecule cell at the Gamma ( Γ ) point, revPBE0-D3 DFTcalculations employing cp2k for a 96 molecule cell at the Γ point, and MLP for a 96 molecule cell at the Γ point. plemented in Crystal17 . Note that the multiple DFTreferences also provide an estimate on the intrinsic er-rors in these DFT calculations due to the choices on the specific hybrid functionals, numerical settings and theuse of different software packages.
Lattice energies and densities . In Fig. 2 we show thecomparison between the lattice energies (left panel) andthe densities (right panel) of the 53 ice phases. Notethat for the lattice energies different theories or exper-iments have different baselines, so for each set of calcu-lation or measurement we use the energies the ice Ih/XIstructure as a reference. In general, we find an excellentagreement between the MLP predictions and the all fourreferences, particularly for the densities. In particular,the differences between the MLP results and the ab ini-tio references are on par with the DFT differences intro-duced by the details of the first principles calculations.For instance, the Pearson correlation coefficient R be-tween the MLP densities and the cp2k values is . and the Root Mean Square Error (RMSE) is . g/mL,and the corresponding metrics between cp2k and vasp are . and . g/mL. For the lattice energies, the R and RMSE between the MLP and the cp2k predictionsare . and meV/H O, respectively, compared with . and meV/H O between the cp2k and the vasp values.
Phonon density of states . The curvature of the poten-tial energy surface around a local minimum relates tothe harmonic frequencies at which the atoms in a crystalvibrate. To investigate the performance of the MLP forthis quantity, we have calculated the phonon frequenciesfor the considered ice structures using the MLP as wellas for a subset of the ice structures using revPBE0-D3DFT calculations with both vasp and cp2k . Figure 3provides a detailed comparison of the phonon densityof states (DOS) for a structure that represent the ice Ihphase and its proton ordered counterpart, XI, showingexcellent agreement in both the low-energy region, cor-responding to long-range dispersive crystal vibrations,and the high-energy region, corresponding to localizedmolecular vibrations. The small shift of low frequencyphonons may be induced by the lack of long-range inter-actions of the MLP. A comparison across all structurescan be found in the Fig. 5 of Methods, and it providesremarkable agreement between the MLP and the firstprinciples methods across all structures in the entireenergy range of vibrations.
III. DISCUSSION
The similarities between the local environments insolid and liquid phases shed light on the structure ofliquid water. There have been many efforts to developa molecular understanding of water, in terms of ori-entational and translational order [4], hydrogen bondnetworks [9] and spontaneously forming dendritic voids[10]. Our approach of using local environments observedin ice as landmark points is a new way of interpretingliquid water as a mixture of ice structures.On the flip side, the similarity also suggests that theliquid and ice structures are distinct in Fig. 1 (a) notbecause of the difference in local environments, but dueto the the presence of long-range order. The conclusionthat liquid water contains all the ice environments ex-plains why the MLP trained on liquid describes the icephases well. This generalization is not specific to thisMLP. Indeed, many water models, such as the coarsegrained mW [23] model, the empirical water modelsSPC [24] and the TIPnP series [25, 26], the polarizableAMOEBA [27], and the MB-pol water potential [28, 29]that are fitted to ab initio reference, qualitatively repro-duce large parts of the phase diagram [30, 31], despitehaving been developed primarily to simulate the liq-uid phase. In particular, MB-pol correctly reproducethe properties of water from the gas to the condensedphases [29].In addition, the general and agnostic comparison ofthe local environments can be easily extended to studyamorphous ice [32], interfaces, and water under con-finement [33]. It is also interesting to investigate hownuclear quantum fluctuations [33–36] influence the dis-tribution of the atomic environments in various phases.Furthermore, our results illustrate the immensepromise of employing MLPs in materials modelling. Forinstance, the MLP used here provides an accurate de-scription of the static and vibrational properties of icephases at a fraction of the cost of the correspondingDFT calculations. For example, the DFT calculationfor the Γ -point phonons of a 4-molecule structure takes CPU hours, and that for a 52-molecule structuretakes , CPU hours, compared with just a few min-utes for both on a laptop using the MLP. Besides, thefact that one can only train a MLP on one liquid phaseand apply the potential to other phases evince the ex-tend of its “extrapolability”, which significantly facili-tates the constructions of the potentials. It is worthnoting that, the MLP for water is stable enough to runMD and PIMD for all the phases. Last but not least, us-ing MLPs as tools for comparing atomic environmentsoffers a new approach of analyzing complex atomic sys-tems in an agnostic and general manner. Our analysisis, of course, not restricted to the chosen DFT level andcan be extended to incorporate new developments inthe field of ab initio methods [37–40].
IV. CONCLUSIONS
To summarize, we compare the local environmentsin various crystalline ice phases and liquid water, usingtwo ML-based approaches. We demonstrate that liq-uid water contains all the local atomic environments indiverse ice phases. Our conclusion provide a new andfundamental perspective on the understanding of liquidwater and ices, and guides future efforts for modelingwater.
V. METHODSA. SOAP representations for atomicenvironments
Numerous representations of atomic environmentshave been developed [14, 41–43], and here we use theSOAP representation [14]. SOAP encodes the local en-vironment X around a central atom using a smoothatomic density function ρ α X ( r ) = (cid:88) | r i | The initial 57 structures are based on an extensivesurvey of ice [13] that generated 15,859 configurations,by exploiting the isomorphism between ice, experimen-tally known zeolites and theoretically-enumerated four-connected SiO networks. The resulting ice-like con-figurations were subsequently locally relaxed, before ageneralized convex hull construction [12] was employedto screen for the ice structures that may be stable undercertain thermodynamic conditions. These include theexperimentally known phases of ice except ice IV, whichwe add back into the selection. Fig. 4 shows the PCAmap of the locations of the selected phases. Notably,many ice phases come in pairs of a low-temperatureproton-ordered form and a higher-temperature proton-disordered form: Ih and XI, III and IX, V and XIII, VIand XV, VII and VIII, and XII and XIV. In this workwe focus on the particular proton-ordered realizationsof these phases made available with Ref. [13]. C. Initial geometry optimization of the icestructures using HSE-3c In Ref. [13] the ice structures were optimized at thePBE DFT level of theory using a coarse k -point gridand plane wave basis, trading accuracy for computa-tional efficiency, For this study, we have therefore per-formed well-converged geometry optimization for thestructures, by running a few cycles of local geometry op-timizations followed by identifying and imposing crys-tal space group symmetries. These local optimizationsare performed with the screened exchange hybrid func-tional HSE-3c [46] using tight optimization thresholdsas implemented in Crystal17 [47, 48]. The Brillouinzone is sampled with a Γ -centered Monkhorst-Pack gridthat has been converged individually for every systemto yield a lattice energy accuracy well below 1 meV.HSE-3c has been shown to yield excellent molecular andintermolecular geometries as well as good noncovalentinteraction energies [49–51], and in particular, suitablefor water and ices. The refined geometries correspondto classical 0 K structures without external pressure andare provided in the Supplemental Information. D. Geometry optimization using VASP The geometries of all ice structures were further re-fined with revPBE0-D3 using the VASP package [52,53]. The equilibrium volumes critically depend on theenergy cutoff, and we used a relatively high value of eV to obtain converged results. The k -point sam-pling grids were the same as those determined for theHSE-3c calculations described above. Structures wereconstrained to their initial symmetries throughout thegeometry optimization, and convergence was achievedwith forces below − eV/Å and stress components be-low − GPa. E. Geometry optimization using cp2k We computed the equilibrium densities and the lat-tice energies of the ice structures using the cp2k code [54] with the revPBE0-D3. The computational de-tails of the calculations are identical with Ref. [11, 55],although the planewave cutoff energy was increased to 800 eV , to obtain smooth volume-energy curves. Weattached the cp2k input file in the Supplemental Infor-mation. Despite a considerable amount of effort, the ge-ometry optimization for 3 structures did not converge toreasonable values, so these calculations were discarded. F. Phonon calculations For the phonon calculations using the MLP, we firstcomputed the Hessian matrix for the 53 geometry-optimized ice phases using finite displacements of 0.01Åof each atom from its equilibrium position along x , y and z axes. Then the Hessian matrix was diagonal-ized to obtain the phonon frequencies as the squareroot of the eigenvalues. We performed those phononcalculations for the ice systems in both their originalcell taken from Ref. [13], and in supercells of this origi-nal cell, obtained by repeating the original cell along allthree crystallographic directions so that each dimensionof the supercell is longer than 8 Å.The phonon calculations using cp2k follow the sameapproach, and the DFT settings are identical to thosein Sec. V E. Presumably due to numerical issues of thespecific DFT setup that we used (e.g. cp2k only sup-ports Γ -point sampling for hybrid functionals), a num-ber of the ice phases contain imaginary phonons at the Γ -point, even after several rounds of geometry optimiza-tion. We discarded the CP2K phonon DOS for thesephases, and only show the ones with real frequencies inFig. 5.The VASP phonon calculations were performed usingthe structures optimized with VASP with the same pa-rameters described above. We used the finite displace-ment method [56] in conjunction with nondiagonal su-percells [57], and commensurate k -point grids were usedto sample the electronic Brillouin zones of the super-cells. The Hessian of a given nondiagonal supercell wascalculated by displacing each atom from its equilibriumposition by . Å in symmetry-inequivalent directionsand calculating the force constants by finite differences.The dynamical matrix for a given q -point grid of the vibrational Brillouin was determined by combining theresults from multiple nondiagonal supercell calculationsas described in Ref. [57]. The resulting dynamical ma-trix was diagonalized to obtain the phonon frequenciesand eigenvectors.In all sets of phonon calculations, imaginary phononsappear in multiple structures at various q -points in theBrillouin zone. This reflects the fact that the protonsin many ice structures are disordered, and when we at-tempt to model them as periodic ordered structures us-ing the unit cells from Ref. [13], we are artificially con-straining them to a saddle point of the potential energysurface rather than to a local minimum. In some of thestructures, instabilities appear even at the Γ -point, andthis is caused by the symmetrization step in preparingthe structures, which again can place them at a saddlepoint. The imaginary phonons in this case break someof the imposed symmetries to lower the energy. Theseproblems can be resolved by replicating the original sim-ulation cell and re-relaxing the atomic positions of thesupercell to allow for the appearance of disorder thatlowers the overall energy, or by re-relaxing the primi-tive cell without imposing symmetry in the case of Γ -point phonons. This additional step is computationallytrivial for the MLP calculations, but computationallyextremely costly for the hybrid functional calculationsusing DFT. As a consequence, it is computationally pro-hibitive to accurately calculate the phonon densities ofstates for all ice structures at the DFT level, and weonly consider a subset in Fig. 5. Data availability The datasets and the Pythonnotebook for analysis are included in the Supplemen-tal Information. The ASAP code is available at: https://github.com/BingqingCheng/ASAP [1] C. G. Salzmann, The Journal of chemical physics ,060901 (2019).[2] E. Brini, C. J. Fennell, M. Fernandez-Serra, B. Hribar-Lee, M. Lukšič, and K. A. Dill, Chemical Reviews ,12385 (2017).[3] D. Marx, M. E. Tuckerman, J. Hutter, and M. Par-rinello, Nature , 601 (1999).[4] J. R. Errington and P. G. Debenedetti, Nature , 318(2001).[5] B. Santra, R. A. DiStasio Jr, F. Martelli, and R. Car,Molecular Physics , 2829 (2015).[6] B. Santra, J. Klimeš, D. Alfè, A. Tkatchenko, B. Slater,A. Michaelides, R. Car, and M. Scheffler, Physical re-view letters , 185701 (2011).[7] M. Fitzner, G. C. Sosso, S. J. Cox, and A. Michaelides,Proceedings of the National Academy of Sciences ,2009 (2019).[8] L. Del Rosso, M. Celli, F. Grazzi, M. Catti, T. C. Hansen, A. D. Fortes, and L. Ulivi, Nature Materials ,1 (2020).[9] I. Ohmine and S. Saito, Accounts of chemical research , 741 (1999).[10] N. Ansari, A. Laio, and A. Hassanali, The Journal ofPhysical Chemistry Letters , 5585 (2019).[11] B. Cheng, E. A. Engel, J. Behler, C. Dellago, andM. Ceriotti, Proceedings of the National Academy ofSciences , 1110 (2019).[12] A. Anelli, E. A. Engel, C. J. Pickard, and M. Ceriotti,Physical Review Materials , 103804 (2018).[13] E. A. Engel, A. Anelli, M. Ceriotti, C. J. Pickard, andR. J. Needs, Nature Communications , 2173 (2018).[14] A. P. Bartók, R. Kondor, and G. Csányi, Phys. Rev.B , 184115 (2013).[15] S. De, A. P. Bartók, G. Csányi, and M. Ceriotti, Phys.Chem. Chem. Phys. , 13754 (2016).[16] J. G. Brandenburg, T. Maas, and S. Grimme, The Figure 5. A comparison between the phonon spectra at the Gamma point for ice phases computed using revPBE0-D3 DFTcalculations employing VASP , revPBE0-D3 DFT calculations employing cp2k , MLP using the original supercell, and MLPusing replicated cell (MLP(L)). For the VASP and cp2k calculations, only those that exhibited no imaginary phonons areshown here. Journal of Chemical Physics , 124104 (2015).[17] Y. Zhang and W. Yang, Phys. Rev. Lett. , 890 (1998).[18] L. Goerigk and S. Grimme, Physical Chemistry Chem-ical Physics , 6670 (2011).[19] S. Grimme, A. Hansen, J. G. Brandenburg, andC. Bannwarth, Chem. Rev. , 5105 (2016).[20] T. E. Markland and M. Ceriotti, Nature Reviews Chem-istry , 1 (2018).[21] J. Behler and M. Parrinello, Physical Review Letters , 146401 (2007).[22] J. Behler, The Journal of chemical physics , 170901(2016).[23] V. Molinero and E. B. Moore, Journal of PhysicalChemistry B , 4008 (2009).[24] H. Berendsen, J. Grigera, and T. Straatsma, Journalof Physical Chemistry , 6269 (1987).[25] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura,R. W. Impey, and M. L. Klein, Journal of ChemicalPhysics , 926 (1983).[26] M. W. Mahoney and W. L. Jorgensen, Journal of Chem-ical Physics , 8910 (2000).[27] M. L. Laury, L.-P. Wang, V. S. Pande, T. Head-Gordon,and J. W. Ponder, The Journal of Physical ChemistryB , 9423 (2015).[28] V. Babin, G. R. Medders, and F. Paesani, Journal ofChemical Theory and Computation , 1599 (2014).[29] S. K. Reddy, S. C. Straight, P. Bajaj, C. Huy Pham,M. Riera, D. R. Moberg, M. A. Morales, C. Knight,A. W. Götz, and F. Paesani, The Journal of chemicalphysics , 194504 (2016).[30] T. Yagasaki, M. Matsumoto, and H. Tanaka, Journalof Physical Chemistry B , 7718 (2018).[31] D. Dhabal, C. Chakravarty, V. Molinero, and H. K.Kashyap, Journal of Chemical Physics , 214502(2016).[32] D. T. Limmer and D. Chandler, Proceedings of the Na-tional Academy of Sciences , 9413 (2014).[33] M. Rossi, M. Ceriotti, and D. E. Manolopoulos, Thejournal of physical chemistry letters , 3001 (2016).[34] D. Marx, M. E. Tuckerman, and M. Parrinello, Journalof Physics: Condensed Matter , A153 (2000).[35] X.-Z. Li, B. Walker, and A. Michaelides, Proceedingsof the National Academy of Sciences , 6369 (2011).[36] F. Paesani, S. Iuchi, and G. A. Voth, The Journal ofchemical physics , 074506 (2007).[37] N. Mardirossian and M. Head-Gordon, J. Chem. Phys. , 214110 (2016).[38] Y. Wang, P. Verma, L. Zhang, Y. Li, Z. Liu, D. G. Truh-lar, and X. He, Proceedings of the National Academyof Sciences , 2294 (2020).[39] M. Riera, E. Lambros, T. T. Nguyen, A. W. Götz, andF. Paesani, Chemical science , 8211 (2019).[40] K. Sharkas, K. Wagle, B. Santra, S. Akter, R. R. Zope,T. Baruah, K. A. Jackson, J. P. Perdew, and J. E. Per-alta, Proceedings of the National Academy of Sciences , 11283 (2020).[41] J. Behler and M. Parrinello, Phys. Rev. Lett. , 146401(2007).[42] F. A. Faber, A. S. Christensen, B. Huang, and O. A.Von Lilienfeld, The Journal of chemical physics ,241717 (2018). [43] A. Sadeghi, S. A. Ghasemi, B. Schaefer, S. Mohr, M. A.Lill, and S. Goedecker, The Journal of chemical physics , 184118 (2013).[44] L. Himanen, M. O. J. Jäger, E. V. Morooka, F. Fed-erici Canova, Y. S. Ranawat, D. Z. Gao, P. Rinke, andA. S. Foster, arXiv e-prints , arXiv:1904.08875 (2019),arXiv:1904.08875 [cond-mat.mtrl-sci].[45] B. Cheng, “Asap: Automatic selection and predictiontools for materials and molecules,” .[46] J. G. Brandenburg, E. Caldeweyher, and S. Grimme,Phys. Chem. Chem. Phys. , 15519 (2016).[47] A. Erba, J. Baima, I. Bush, R. Orlando, and R. Dovesi,J. Chem. Theory Comput. , 5019 (2017).[48] R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, L. Maschio, M. Rérat, S. Casassa,J. Baima, S. Salustro, and B. Kirtman, WIREs Com-put. Mol. Sci. , e1360 (2018).[49] E. Caldeweyher and J. G. Brandenburg, J. Phys.: Con-dens. Matter , 213001 (2018).[50] S. Rösel, H. Quanz, C. Logemann, J. Becker, E. Mossou,L. Canadillas-Delgado, E. Caldeweyher, S. Grimme,and P. R. Schreiner, J. Am. Chem. Soc. , 7428(2017).[51] L. Doná, J. G. Brandenburg, and B. Civalleri, J. Chem.Phys. , 121101 (2019).[52] G. Kresse and J. Furthmüller, Phys. Rev. B , 11169(1996).[53] G. Kresse, Phys. Rev. B , 169 (1996).[54] G. Lippert, J. Hutter, and M. Parrinello, TheoreticalChemistry Accounts , 124 (1999).[55] O. Marsalek and T. E. Markland, J. Phys. Chem. Lett. , 1545 (2017).[56] K. Kunc and R. M. Martin, Phys. Rev. Lett. , 406(1982).[57] J. H. Lloyd-Williams and B. Monserrat, Phys. Rev. B92