ARTICLE DOI: 10.1038/s41467-018-04618-6 OPEN Mapping uncharted territory in ice from zeolite networks to ice structures 1 2 2 3,4 1 Edgar A. Engel , Andrea Anelli , Michele Ceriotti , Chris J. Pickard & Richard J. Needs Ice is one of the most extensively studied condensed matter systems. Yet, both experi- mentally and theoretically several new phases have been discovered over the last years. Here we report a large-scale density-functional-theory study of the conﬁguration space of water ice. We geometry optimise 74,963 ice structures, which are selected and constructed from over ﬁve million tetrahedral networks listed in the databases of Treacy, Deem, and the International Zeolite Association. All prior knowledge of ice is set aside and we introduce “generalised convex hulls” to identify conﬁgurations stabilised by appropriate thermodynamic constraints. We thereby rediscover all known phases (I–XVII, i, 0 and the quartz phase) except the metastable ice IV. Crucially, we also ﬁnd promising candidates for ices XVIII through LI. Using the “sketch-map” dimensionality-reduction algorithm we construct an a priori, navigable map of conﬁguration space, which reproduces similarity relations between structures and highlights the novel candidates. By relating the known phases to the tractably small, yet structurally diverse set of synthesisable candidate structures, we provide an excellent starting point for identifying formation pathways. 1 2 TCM Group, Cavendish Laboratory, J J Thomson Avenue, Cambridge CB3 0HE, UK. Laboratory of Computational Science and Modeling, Institute of Materials, École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland. Department of Materials Science and Metallurgy, University of Cambridge, 27 Charles Babbage Road, Cambridge CB3 0FS, UK. Advanced Institute for Materials Research, Tohoku University, 2-1-1 Katahira, Aoba, Sendai 980-8577, Japan. Correspondence and requests for materials should be addressed to E.A.E. (email: email@example.com) NATURE COMMUNICATIONS (2018) 9:2173 DOI: 10.1038/s41467-018-04618-6 www.nature.com/naturecommunications 1 | | | 1234567890():,; ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04618-6 ce is a complex system of interest across much of science, 10 500 ranging from astrophysics to biology. On the Earth’s surface Iand in its atmosphere, it plays a central role in determining 9 climate and in countless natural processes and technological applications. Ice is also a key constituent of the Earth’s crust and mantle. Its phase diagram and properties have been investigated across a wide range of temperatures and pressures by experi- mentalists and theoreticians alike. A total of 18 crystalline ice phases have been formed under 200 1 2 various conditions , 7 of which are metastable . In addition, a 3–14 number of hypothetical ice phases have been predicted and characterised using computer simulations. All of these phases are molecular crystals that fulﬁl the “Bernal–Fowler ice rules” and form four-connected networks. In most ice phases, the distinct 3 0 ways of dressing the oxygen sublattice with hydrogen atoms 34 56 789 10 within the ice rules (the so-called “proton-orderings”) are quasi- SiO energetically degenerate . Theoretical studies have also suggested structures of water ice under ultrahigh pressures of up to many Fig. 1 Isomorphism between SiO and H O polymorphs. Correlation 2 2 17,18 terapascals and its eventual decomposition . between the average ring sizes, r, of SiO and H O polymorphs. More than 2 2 The phase diagram of ice has recently received renewed 1/3 of the ice polymorphs retain the ring statistics of their counterpart SiO interest, ﬁrst, because the theoretical discovery of the s-III network 13 19 clathrate hydrate and the experimental description of ice XVII 14,20,21 and of two-dimensional forms of ice have demonstrated Results that our understanding of ice is far from complete. Second, it has Exploring conﬁguration space. The strong isomorphism become apparent that the nucleation and melting of ice are between ice and silica networks has previously been explored in 22,23 complex processes in which metastable ice phases play a role . ref. and arises because both silica and water preferentially form In classical nucleation theory, an interfacial free energy advantage four-connected networks composed of corner-sharing tetrahedral of a few percent will lead to preferential nucleation of metastable units. The basic building blocks of silica and water ice are so phases with free energies up to around 10 meV/H O above the similar that it is even possible to form silica/water hetero- stable phase . networks in which silicate oligomers form part of the hydrate Despite valiant efforts using structure searching methods such 35,36 lattice . as ab initio random structure searching (AIRSS) , to our There is a vast literature on four-connected structures, knowledge no comprehensive study of (meta)stable ice phases including an atlas describing the underlying networks of porous and their formation has been published to date. The problem is crystalline zeolites and a number of very large databases of two-fold: ﬁrst, the enormous conﬁguration space must be theoretically enumerated networks, such as the databases of 38 39 explored in a reasonably comprehensive manner. Second, in Treacy and Deem . Graph network enumeration has pre- order to render the structure search relevant to experiment, the viously been applied to crystal structure prediction and, in large number of theoretical (meta)stable structures generated in 2 3 41,42 particular, to sp - and sp -carbon . The above databases have the process must be reduced to those that can be formed proven to be a valuable resource in searching for sp allotropes of experimentally. This reﬁnement must be a priori and quantitative. 43 carbon and constitute a comprehensive source of four- Finally, different stabilising factors—such as the absorption of connected networks from which topologically distinct phases of guest molecules —can be investigated further, and methods such ice can be constructed and geometry optimised to the respective 27–30 as forward ﬂux sampling and enhanced sampling metady- associated local minimum energy structures using conjugate 31–33 namics may be used to identify possible synthetic pathways. gradient methods. Recently, the search for computationally stable 44 37 This work aims for a comprehensive study of crystalline ice ultralow-density ices on the basis of the atlas of zeolites has phases, focussing on the exploration of conﬁguration and the hinted at the potential of this approach, despite its more limited reduction of the resulting intractably large amount of structure scope and despite only considering stabilisation under (effective data to a small number of structures, which are likely to be negative) pressure. The search for (meta)stable ice phases is accessible experimentally. In the ﬁrst section of the Results we facilitated by the strong correspondence between zeolite and ice exploit the isomorphism between ice and silica networks to structures. Figure 1 shows the strong correlation between the explore the relevant parts of the conﬁguration space of ice using average ring sizes of SiO structures and their counterpart H O 2 2 databases of theoretically enumerated, four-connected networks. polymorphs after geometry optimisation, indicating that structu- In the second section of the Results we rationalise the resultant rally distinct SiO networks generally translate into structurally structural data on the basis of purely energetic considerations and distinct H O networks. thereby identify structures that can be stabilised under pressure. The large size of the databases of hypothetical zeolites By design, this approach cannot identify structures stabilised by necessitates some preselection of structures. Tribello et al. show thermodynamic and kinetic constraints other than pressure, such that the energies and densities of low-density SiO networks and as temperature, electric ﬁelds, concentrations of guest molecules, their counterpart H O networks are correlated, but this correla- etc. In the third section of the Results we overcome this limitation tion does not carry across to structures of densities comparable to by developing a “generalised convex hull” (GCH) construction. and higher than that of ice Ih (see Supplementary Note 1 and Moreover, we use the sketch-map algorithm to construct a Fig. 1). Consequently, neither SiO lattice energies nor densities navigable map of the conﬁguration space of ice. This primarily can be used for preselection. Since all known ice phases (with the serves as an aide in developing an intuitive understanding of exception of ice V/XIII) have unit cells containing no more than structural relationships. However, it also shows potential for 16 molecules, applying a cutoff to the unit cell size provides a helping to identify formation pathways for new candidate ice reasonable method for preselection, which can be improved phases. 2 NATURE COMMUNICATIONS (2018) 9:2173 DOI: 10.1038/s41467-018-04618-6 www.nature.com/naturecommunications | | | H O # of structures NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04618-6 ARTICLE systematically by including structures with larger unit cells. In Static practice, we preselect only networks with unit cell volumes of no XIII Harmonic more than 800 Å and without 3-rings, which would normally II induce excessive strain in an ice structure. Out of the 331,172 (Deem) and 5,389,408 (Treacy) zeolites, this leaves 74,731 struc- tures. This selection contains duplicates since the databases are Q III not mutually exclusive. Low-density structures with low SiO IX lattice energies are added back in by including the experimentally synthesised zeolite networks from the IZA database (see also 15 Supplementary Fig. 2). Geometry optimisation of the resulting 74,963 structures using ﬁrst-principles quantum-mechanical methods is viable. However, HS-III at this stage only rough lattice energies are required to identify the XVII low-energy sectors of conﬁguration space. The deﬁnition of low energies is provided by the differences between the lattice energies 0.80 0.90 1.00 1.10 1.20 1.30 of different proton-orderings and between the quantum vibra- tional corrections of different structures, which are both of ρ (g/cm ) the order of 10 meV/H O . Benchmarking against more accurate Fig. 2 Energy-density convex hull. PBE-DFT static lattice energies (red) and density-functional-theory results using the PBE exchange- free energies including harmonic vibrations (blue) relative to ice Ih for correlation functional (PBE-DFT) (see Supplementary Note 5), known ice phases (blue labels) and energetically competitive phases (black which are further benchmarked against results obtained using the 48 labels). The labels of the novel energetically competitive phases correspond rPW86-vdW2 functional in Supplementary Note 6, shows that 49 to the numbering scheme in Fig. 3. The energy-density convex hulls at the energies from ReaxFF force ﬁelds are sufﬁcient for this purpose static lattice and harmonic vibrational levels are indicated by red and blue (see Supplementary Fig. 5). After removing high-energy conﬁg- solid lines, respectively urations, the geometries of the remaining structures are reﬁned using PBE-DFT. Removing duplicates leaves 15,869 distinct Ih/XI, II, III/IX, V/XIII, VII/VIII and X phases of ice. Moreover, it structures. identiﬁes the structure that has since been identiﬁed experimen- 19,50 tally as the porous ice XVII . This clearly demonstrates the potential of our structure searching approach. However, not all Phase stability and characterisation of structures. The large known ice phases are classiﬁed as synthesisable, which highlights pool of candidate structures highlights the central challenge of the limitations of the established convex-hull approach: it fails to computational structure searches: the number of theoretical identify synthesisable metastable structures (such as ice IV and (meta)stable conﬁgurations that can be constructed increases XII/XIV) and structures that can be stabilised and made exponentially with system size, but only those that can be synthesisable by thermodynamic and kinetic constraints other observed experimentally are of interest. Their selection must take than pressure (such as XVI that initially forms by absorption of into consideration the uncertainty in the computational frame- H guest molecules). These limitations will be addressed in the work, the possibility of kinetic and/or surface effects promoting 2 following. the formation of metastable phases and the (de)stabilisation of In addition, the ice counterparts of the zeolites with network phases by different thermodynamic boundary conditions, such as codes IRR, IWV, SGT and DDR and three hypothetical zeolites pressure. are identiﬁed as prime candidates for stabilisation by varying the To identify the polymorphs that are most likely to form at system density. The counterparts of the IRR and IWV (not shown different pressures, we ﬁrst consider a well-established approach in Fig. 2), 207_1_4435 and DDR zeolites (labelled 1 and 17 in based on a convex-hull construction. The convex hull of energy Fig. 2) are excellent candidates for stabilisation under negative (as a proxy for free energy) as a function of density, E (ρ), is ch pressure or by inclusion of guest molecules. The DDR counter- formed by structures that are stable against decomposition into part, in particular, has previously been proposed as a possible two or more structures with lower average energy at the same clathrate hydrate . The IRR and IWV counterparts exhibit average density and the so called “tie lines” that connect them. In substantially lower densities than the known CS-I, CS-II and HS- the absence of kinetic effects, the only phases that can be observed 51–54 III clathrate hydrates , suggesting that they may only become by manipulating the density of the system (for example through stable at large negative pressures. Conversely, the counterparts of pressure) are exactly those that constitute the vertices of the the PCOD8172143 and 11_2_15848 zeolites (labelled 15 and 18 energy-density convex hull (see Fig. 2). In analogy with the in Fig. 2) may be stabilised under positive pressure (also see Bell–Evans–Polanyi principle, which states that highly exothermic Supplementary Fig. 4). chemical reactions have low activation energies, the stability of a When comparing structures whose stabilities lie within a few given metastable structure can be assessed by the free energy of meV/H O of each other, anharmonic quantum nuclear effects decomposition into stable structures. We refer to this as the 2 (QNE) must be accounted for, as highlighted by the stabilisation “dressed energy”. Plainly put, the proximity of a metastable of ice Ih with respect to Ic by anharmonic quantum nuclear structure to the convex hull is a measure of its stability. The vibrations , as well as by effects of similar magnitude observed in “dressed energy” is calculated by subtracting the convex-hull other H-bonded crystals . Anharmonic QNE in particular energy at the corresponding density ρ from the lattice energy E stabilise ice XVII and the HS-III clathrate by a few meV/H O (as a proxy for free energy), E = E − E (ρ). Based on E ,we 2 dr ch dr with respect to Ih (see Table SIII). Their resultant zero pressure reﬁne the selection of ice structures as speciﬁed in Computational free energies exceed that of Ih by only 6.8 and 7.8 meV/H Oat methods. Ultimately, only structures with E less than 10 meV dr the PBE-DFT level, respectively. The relative stability of the are retained, for which kinetic, entropic and/or surface effects counterparts of the PCOD8172143 and 11_2_15848 zeolites with may plausibly lead to preferential formation during nucleation. respect to Ih, on the other hand, is affected very little (see Setting aside all prior knowledge of ice, this procedure Supplementary Note 4 for further detail). identiﬁes the theoretical i, 0 and quartz phases and the known NATURE COMMUNICATIONS (2018) 9:2173 DOI: 10.1038/s41467-018-04618-6 www.nature.com/naturecommunications 3 | | | E (meV/H O) 2 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04618-6 12 30 i Ih / XI III/IX Ic V / XIII XII / XIV VI / XV XVII 31 Square 25 II CS-II 1 1.5 2 19 g/cm HS-III VII / VIII GCH1 21 0 0.2 0.4 0.6 0.8 GCH3 ΔE [eV] dr Fig. 3 Sketch map of the structural similarity of 15,869 distinct PBE-DFT geometry-optimised ice structures. The sketch-map coordinates correlate strongly with density and conﬁgurational energy but ultimately measure abstract structural features, which leaves their numerical value without intuitive meaning. They are therefore not shown. Instead the density and static lattice energy of each structure is encoded by the size and colour of the respective point on the map. Known ice phases are labelled in blue. The 34 new candidates are labelled in black and numbered in order of increasing dressed energy relative to the GCH3. Their atomic structures are shown to highlight their structural diversity Using machine-learning to navigate the structural landscape. kernel-induced distance, we apply the sketch-map algorithm to An analysis based on the energy-density convex hull as in the obtain a two-dimensional representation that reproduces as second section of the Results identiﬁes candidate structures that accurately as possible the (non-linearly transformed) distance can be stabilised by pressure. However, this does not address between each pair of structures. The construction and its several crucial issues: (a) obtaining a global picture of conﬁg- parameters were designed to assess the oxygen lattice while being uration space from which one can gather an intuitive under- insensitive to proton-disorder and hydrogen-bonding defects. standing of the relations between different polymorphs; (b) The resulting map is shown in Fig. 3 and provides a much- assessing the effectiveness of the structure search, identifying needed global picture of the lie of the land. Notably, it is spanned more or less obvious “gaps”; and (c) selecting candidates stabi- by collective coordinates measuring abstract structural features, lised by thermodynamic constraints other than pressure, such as which in general cannot be related to single conventional absorption of guest molecules, electric ﬁelds, etc. All of these observables such as density in a meaningful way. Consequently, problems can be tackled effectively within a framework that their numerical values are not shown in Fig. 3. borrows ideas from the machine-learning community. Points (a) Several observations highlight the heuristic value of such a and (b) are addressed by constructing an abstract, unbiased and representation: (1) The positions on the map correlate well with general two-dimensional representation of conﬁguration space in both density and lattice energy (see Supplementary Note 2 and terms of the similarity relations between structures. Point (c) is Supplementary Fig. 5); (2) Structures related by proton-disorder, addressed by generalising the conventional convex hull such as Ih/XI, III/IX and VII/VIII, are clustered together; (3) construction. Structures related by stacking disorder, such as Ih, Ic and Isd, are The ﬁrst key ingredient of an intuitive representation of clustered together; (4) The spread in energy at a given point on conﬁguration space is a measure of the similarity of different the map is comparable to the energy scale of stacking defects and conﬁgurations. We use the smooth overlap of atomic positions H-bonding defects. H-bonding defects and different proton- (SOAP) kernel combined with an entropy-regularised matching orderings develop during the geometry optimisation of the ice (REMatch) approach . This captures the fundamental symme- structures, which (in analogy with their SiO parent structures) tries of the problem, such as invariance to alternative representa- are initialised with bond-centred protons. tions of the same periodic structure, particle labelling and rigid Furthermore, the general structure of the map is consistent rotations and translations of the atomic coordinates. Based on the with the strategy we followed to construct our set of structures. 4 NATURE COMMUNICATIONS (2018) 9:2173 DOI: 10.1038/s41467-018-04618-6 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04618-6 ARTICLE The upper portion of the map, corresponding to tetrahedral ices Furthermore, we identify 34 new conﬁgurations that are excellent and silica-like networks is densely sampled, with structures candidates for experimental formation and that we propose as clustered in partially overlapping regions. The lower part of the candidates for ices XVIII through LI. Among them are, in map corresponding to very dense (e.g., ice X) and very open particular, the ice counterparts of the DDR, SGT and NON structures (e.g., those originating from the IZA zeolite data set) is zeolites, which were previously suggested as promising candidates sparse. At high density, this sparsity results from the increasing for clathrate hydrates by Tribello et al. , and two structures importance of geometric constraints, which limit structural reminiscent of a high pressure structure with Pbcm symmetry diversity and prohibit the formation of energetically feasible proposed by Hermann et al. . Notably, while the most promising “mixed phases” containing structural patterns from two or more candidates for experimental formation (as indicated by their low-energy conﬁgurations. At low-density, our preselection ordering in Fig. 3) are low-density ice counterparts of different strategy leads to sparse sampling. Sketch map therefore provides zeolite networks, the counterpart of the ITT network, which was indications of the quality of conﬁguration-space sampling, which suggested as the most stable “aeroice” structure below around can be used to focus the structure search on the regions that need −0.4 GPa in ref. , is dynamically unstable at the employed level it most. of theory. For reference, using the rPW86-vdW2 exchange- Finally, structures with low E are projected onto the correlation functional ITT ice is still much less stable than IRR ice dr periphery of the map (see Supplementary Fig. 4c), whereas the proposed as stabilisable in this work (see Supplementary Fig. 1b). central region is largely populated by defective, “mixed phase” It is worth noting that the counterpart of the LTA zeolite structures that lie far from the energy-density convex hull. This (structure 4 in Fig. 3) has also most recently received attention as suggests that a machine-learning-inspired analysis of structures an “ultralow” density clathrate ice in ref. . While the GCH may be used to establish a GCH, which identiﬁes conﬁgurations generally depends on the kernel, our choice of kernel representa- that can be stabilised (and made “synthesisable” ) by the tion is very general and rather unbiased, which is reﬂected by the application of appropriate thermodynamic constraints. weak dependence of this selection of structures on the choice of We deﬁne this GCH construction in analogy with the hyperparameters for the SOAP-REMatch kernel. Notably, the conventional energy-density convex hull in the second section GCH is also remarkably insensitive to the details of the of the Results, but instead of considering E as a function of ρ,we underlying (free) energy calculations. As shown in Supplementary consider E as a function of n variables measuring abstract Note 6, 37 out of the 38 structures are still identiﬁed as GCH structural features, E(ϕ , …, ϕ ). The simplices of the GCH thus vertices when the lattice energies of the structures highlighted in 1 n correspond to structures that are stable with respect to Fig. 3 are computed using the dispersion-corrected rPW86-vdW2 decomposition at a given set of these abstract structural features, exchange-correlation functional instead of the PBE functional, rather than at a particular density. The conventional energy- despite signiﬁcant differences with respect to the PBE lattice density convex hull allowed the identiﬁcation of structures, which energies. The rPW86-vdW2 functional has been shown to be can be stabilised by pressure, because pressure allows the particularly accurate for the known phases of ice . In contrast, manipulation of the density of an ice sample. Conversely, the the energy-density CH depends more sensitively on the choice of GCH allows us to identify structures, which can be stabilised by exchange-correlation functional. imposing thermodynamic and/or kinetic constraints that couple to the abstract structural features. In analogy with E , one can dr ðnÞ then deﬁne a generalised dressed energy E that quantiﬁes the Discussion dr stability of a given conﬁguration subject to constraints that couple The success of the GCH construction in discovering the known to the n structural features ϕ … ϕ . ice phases and clathrate hydrates entirely a priori highlights that, 1 n In practice, a kernel principal-component analysis (KPCA) is although kinetic factors play an important role in determining performed to extract the KPCA descriptors ϕ , which encode the which ice phases are formed in practice, structural and simple key structural features and can be sorted in order of decreasing energetic considerations can provide a great deal of physical importance. Crucially, the KPCA components (unlike the insight. More importantly, it demonstrates that the GCH collective variables deﬁning the highly non-linear sketch-map approach does not simply discern structurally diverse conﬁg- projection) form a vector space in which the notion of convexity urations, but very effectively selects conﬁgurations that can be is well deﬁned. formed in experiment. It thereby provides strong support for the The GCH construction provides a powerful tool for discovery, 34 proposed, new, structurally diverse candidates for ices since one can select conﬁgurations that are both low in energy XVIII–LI. This should spur experimental efforts to ratify our and “extremal” in the sense of structural features described by one predictions. or more KPCA descriptors. By increasing the number n of At this stage, the candidates are embedded in a human- features considered, the screening becomes progressively more readable sketch map of conﬁguration space mainly as an aide in inclusive, since multiple axes of structural diversity are considered developing an intuitive understanding of the relation of the simultaneously. In practice, the selection was further reﬁned by proposed candidates to the known ice phases. However, high- (automatically) eliminating structurally related conﬁgurations, as lighting the phase transitions between the known ice phases discussed in Supplementary Note 2. suggests that proximity on the sketch map is a good indicator for Including three KPCA descriptors in the GCH construction, the existence of a viable transition pathway. In conjunction with we identify 50 structures within 20 meV of the GCH (see Fig. 3), the GCH construction, the sketch map therefore provides a which include all of the known ice phases except ice IV. Ice IV is tractably small and yet structurally diverse set of synthesisable not classiﬁed as synthesisable due to its particularly high lattice candidate structures and a means of identifying end points and energy, which is consistent with the experimental observation suitable reaction coordinates for further investigation of forma- that ice IV is metastable and only forms occasionally upon slow tion pathways, for example, using transition state sampling , 64 27 heating of high-density amorphous ice before annealing to ice III, umbrella sampling , forward ﬂux sampling or enhanced V or VI. sampling metadynamics approaches , which have already pro- 28–32 The 50 structures also include the theoretical i, 0, quartz and ven successful in simulating the nucleation of ice . The square phases, the CS-II clathrate hydrate (which is identical in relation of the KPCA descriptors in the GCH construction to structure to ice XVI) and the HS-III clathrate hydrate. conventional quantities, such as density, vibrational spectra and NATURE COMMUNICATIONS (2018) 9:2173 DOI: 10.1038/s41467-018-04618-6 www.nature.com/naturecommunications 5 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04618-6 concentrations of different types of guest molecules, promises to Data availability. The data that support the ﬁndings of this study are available at the following https://doi.org/10.24435/materialscloud:2018.0010/v1. provide more direct guidance in identifying experimental for- mation pathways. However, this goes beyond the scope of this Received: 22 October 2017 Accepted: 11 May 2018 study. In addition, the approach demonstrated in this work sheds light on the energetics of proton-order/disorder and stacking- disorder, as well as H-bonding and planar defects, and also provides a glimpse of the preferred (quasi-) two-dimensional forms of ice. References In its current state, the biggest limitation of our structure 1. Hobbs, P. V. Ice Physics (Oxford University Press, Oxford, 2010). 2. Petrenko, V. F. & Whitworth, R. W. Physics of Ice (Oxford University Press, search is the preselection cutoff on system size. Relaxing this Oxford, 1999). cutoff will drive the structure search towards completeness, which 3. Singer, S. J. et al. Hydrogen-bond topology and the ice VII/VIII and Ice Ih/XI is the obvious next step. More generally, the connection between proton-ordering phase transitions. Phys. Rev. Lett. 94, 135701 (2005). structural patterns and conﬁgurational energy exposed by the 4. Knight, C. & Singer, S. J. Prediction of a phase transition to a hydrogen bond sketch-map dimensionality reduction suggests an expedient ordered form of ice VI. J. Phys. Chem. B 109, 21040 (2005). 5. Kuo, J.-L. The low-temperature proton-ordered phases of ice predicted by ab recipe for even more extensive database-driven searches. These initio methods. Phys. Chem. Chem. Phys. 7, 3733 (2005). need not be limited to crystalline water ice but could range from 6. Kuo, J.-L. & Kuhs, W. F. A ﬁrst principles study on the structure of ice-VI: other tetrahedrally coordinated systems, such as silica or the static distortion, molecular geometry, and proton ordering. J. Phys. Chem. B carbon allotropes, to liquid, disordered and glassy systems. 110, 3697 (2006). 7. Knight, C. & Singer, S. J. Hydrogen bond ordering in ice V and the transition to ice XIII. J. Chem. Phys. 129, 164513 (2008). 8. Tribello, G. A., Slater, B. & Salzmann, C. G. A blind structure prediction of ice Methods XIV. J. Am. Chem. Soc. 128, 12594 (2006). Geometry optimisations. Ice structures were initially geometry optimised using 9. Russo, J., Romano, F. & Tanaka, H. New metastable form of ice and its role in 49 66 ReaxFF force ﬁelds as parametrised by Raymand et al. and implemented in the the homogeneous crystallization of water. Nat. Mater. 13, 733 (2014). 67 −4 Gulp package until the energies and forces were converged to within 10 eV and 10. Fennell, C. J. & Gezelter, J. D. Computational free energy studies of a new ice −3 10 eV/Å, respectively. This choice is motivated and justiﬁed in Supplementary polymorph which exhibits greater stability than ice Ih. J. Chem. Theory Note 5. Comput. 1, 662 (2005). First-principles quantum-mechanical geometry optimisations were performed 47 68 11. Svishchev, I. M. & Kusalik, P. G. Quartzlike polymorph of ice. Phys. Rev. B 53, using semi-local PBE-DFT as implemented in the Castep package . The choice R8815 (1996). of density functional is discussed in detail in Supplementary Note 6. The initial 12. Tribello, G. A., Slater, B., Zwijnenburg, M. A. & Bell, R. G. Isomorphism PBE-DFT geometry optimisations were performed with an plane-wave energy between ice and silica. Phys. Chem. Chem. Phys. 12, 8597 (2010). cutoff of 490 eV, Monkhorst-Pack k-point grids of maximum spacing of 2π × 0.07 −1 13. Huang, Y. et al. A new phase diagram of water under negative pressure: the Å , and on-the-ﬂy generated ultrasoft pseudopotentials. The structures in the rise of the lowest-density clathrate s-III. Sci. Adv. 2, e1501010 (2016). second section of the Results were further reﬁned using PBE-DFT calculations with 14. Ji, C., Schusteritsch, G., Pickard, C. J., Salzmann, C. G. & Michaelides, A. Two norm-conserving pseudo-potentials and, ﬁrst, a plane-wave energy cutoff of 490 eV −1 dimensional ice from ﬁrst principles: structures and phase transitions. Phys. and Monkhorst-Pack k-point grids of maximum spacing of 2π × 0.07 Å , then a Rev. Lett. 116, 025501 (2016). cutoff of 800 eV and Monkhorst-Pack k-point grids of maximum spacing of 2π × −1 0.04 Å , and ﬁnally a cutoff of 1200 eV and Monkhorst-Pack k-point grids of 15. Bernal, J. D. & Fowler, R. H. A theory of water and ionic solution, with −1 maximum spacing of 2π × 0.04 Å . The resulting energy differences between particular reference to hydrogen and hydroxyl ions. J. Chem. Phys. 1, 515 frozen-phonon conﬁgurations, atomic positions and residual forces were converged (1933). −4 −5 −4 to within 10 eV/H O, 10 Å and 10 eV/Å, respectively (also see 16. Pauling, L. The structure and entropy of ice and of other crystals with some Supplementary Note 3). randomness of atomic arrangement. J. Am. Chem. Soc. 57, 2680 (1935). Harmonic vibrational modes and frequencies were calculated using a ﬁnite 17. Hama, J. & Suito, K. Physics and Chemistry of Ice (Hokkaido University Press, displacement method. Anharmonic vibrations are calculated using the vibrational Sapporo, 1992). self-consistent ﬁeld approach described in ref. . The 3N-dimensional BO energy 18. Pickard, C. J., Martinez-Canales, M. & Needs, R. J. Decomposition and surface (where N is the number of atoms in the simulation cell) was described by terapascal phases of water ice. Phys. Rev. Lett. 110, 245701 (2013). mapping one-dimensional (1D) subspaces along the harmonic normal mode axes 19. del Rosso, L., Celli, M. & Ulivi, L. New porous water ice metastable at up to large amplitudes of four times the harmonic root-mean-square atmospheric pressure obtained by emptying a hydrogen-ﬁlled ice. Nat. displacements, where anharmonicity is important. The 3N-dimensional BO surface Commun. 7, 13394 (2016). was then reconstructed from the 1D subspaces. The 1D energy surfaces were ﬁtted 20. Algara-Siller, G. et al. Square ice in graphene nanocapillaries. Nature 519, 443 using cubic splines The anharmonic vibrational Schrödinger equation was solved (2015). expanding the vibrational wave function in terms of simple harmonic oscillator 21. Ji, C., Schusteritsch, G., Pickard, C. J., Salzmann, C. G. & Michaelides, A. eigenstates. The inclusion of 25 states for each vibrational degree of freedom was Double-layer ice from ﬁrst principles. Phys. Rev. B 95, 094121 (2017). found sufﬁcient to obtain converged results. 22. Haji-Akbari, A. & Debenedetti, P. G. Direct calculation of ice homogeneous Duplicates were identiﬁed by applying the “crysim” tool from the AIRSS nucleation rate for a molecular model of water. Proc. Natl. Acad. Sci. USA 112, method to the oxygen sublattices. 10582 (2015). 23. Quigley, D. Communication: Thermodynamics of stacking disorder in ice nuclei. J. Chem. Phys. 141, 121101 (2014). 24. Quigley, D., Alfé, D. & Slater, B. Communication: On the stability of ice 0, ice Machine-learning analysis of structural relations. To assess the structural similarity between the conﬁgurations in the database under study, we used a i, and Ih. J. Chem. Phys. 141, 161102 (2014). 25. Pickard, C. J. & Needs, R. J. Ab initio random structure searching. J. Phys. REMatch-SOAP kernel, as implemented in the glosim.py package (http:// cosmo-epﬂ.github.io), with the following choice of hyperparameters controlling the Condens. Matter 23, 053201 (2011). description of atomic environments: 26. Tribello, G. A. & Slater, B. A theoretical examination of known and /src/glosim/glosim.py -n 9 -l 6 -c 5 -g 0.5 hypothetical clathrate hydrate materials. J. Chem. Phys. 131, 024703 (2009). periodic nocenter 1 kernel rematch gamma 0.01 27. Allen, R. J., Frenkel, D. & ten Wolde, P. R. Forward ﬂux sampling-type nonorm schemes for simulating rare events: efﬁciency analysis. J. Chem. Phys. 124, Hydrogen atoms were included in the deﬁnition of the atom-density overlap 024102 (2006). but were not considered as environment centres, so as to de-emphasise proton- 28. Li, T., Donadio, D. & Galli, G. Ice nucleation at the nanoscale probes no man's (dis)order in the deﬁnition of structural similarity. The choice of cutoff radius was land of water. Nat. Commun. 4, 1887 (2013). tuned to achieve a clear separation between the known phases of ice in the 29. Haji-Akbari, A. & Debenedetti, P. G. Direct calculation of ice homogeneous database. nucleation rate for a molecular model of water. PNAS 112, 10582 (2015). The non-linear sketch-map dimensionality reduction scheme (http:// 30. Bi, Y., Cao, B. & Li, T. Enhanced heterogeneous ice nucleation by special sketchmap.org) was then applied to the SOAP kernel measure of similarity for 400 surface geometry. Nat. Commun. 8, 15372 (2017). farthest-point-sampled landmark structures following the procedure described in 31. Quigley, D. & Rodger, P. M. A metadynamics-based approach to sampling ref. and using the following parameters: σ = 0.12, A = 2, B = 4, a = 2, and b = 2. crystallisation events. Mol. Simul. 35, 613 (2009). 6 NATURE COMMUNICATIONS (2018) 9:2173 DOI: 10.1038/s41467-018-04618-6 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04618-6 ARTICLE 32. Giberti, F., Salvalaglio, M. & Parrinello, M. Metadynamics studies of crystal 62. Santra, B. et al. On the accuracy of van der Waals inclusive density-functional nucleation. IUCrJ 2, 256 (2015). theory exchange-correlation functionals for ice at ambient and high pressures. 33. Pipolo, S. et al. Navigating at will on the water phase diagram. Phys. Rev. Lett. J. Chem. Phys. 139, 154702 (2013). 119, 245701 (2017). 63. Dellago, C., Bolhuis, P. & Geissler, P. L. Transition path sampling. Adv. Chem. 34. Ceriotti, M., Tribello, G. A. & Parrinello, M. Simplifying the representation of Phys. 123, 1 (2002). complex free-energy landscapes using sketch-map. Proc. Natl. Acad. Sci. USA 64. Kumar, S., Rosenberg, J. M., Bouzida, D., Swendsen, R. H. & Kollman, P. A. 108, 13023 (2011). The weighted histogram analysis method for free-energy calculations on 35. Emmer, J. & Wiebcke, M. Heteronetwork clathrates with three-dimensional biomolecules. I. The method. J. Comp. Chem. 13, 1011 (1992). mixed silicate-water host frameworks and channel systems. J. Chem. Soc. 65. Laio, A. & Parrinello, M. Escaping free energy minima. Proc. Natl. Acad. Sci. Chem. Commun. 2079–2080 (1994). USA 99, 12562 (2002). 36. Wiebcke, M. Structural links between zeolite-type and clathrate hydrate-type 66. Raymand, D. et al. Water adsorption on stepped ZnO surfaces from MD materials. J. Chem. Soc. Chem. Commun. 1507–1508 (1991). simulation. Surf. Sci. 604, 741 (2010). 37. Baerlocher, C., Meier, W. M. & Olson, D. H. Atlas of Zeolite Framework Types 67. Gale, J. D. & Rohl, A. L. The General Utility Lattice Program (GULP). Mol. (Elsevier, Amsterdam, 2007). Simul. 29, 291 (2003). 38. Treacy, M. M. J., Rivin, I., Balkovsky, E., Randall, K. H. & Foster, M. D. 68. Clark, S. J. et al. First principles methods using CASTEP. Z. Kristallogr. 220, Enumeration of periodic tetrahedral frameworks. II. Polynodal graphs. 567 (2005). Microporous Mesoporous Mater. 74, 121 (2004). 69. Monserrat, B., Drummond, N. D. & Needs, R. J. Anharmonic vibrational 39. Earl, D. J. & Deem, M. W. Toward a database of hypothetical zeolite properties in periodic systems: energy, electron-phonon coupling, and stress. structures. Ind. Eng. Chem. Res. 45, 5449 (2006). Phys. Rev. B 87, 144302 (2013). 40. Winkler, B., Pickard, C. J., Milman, V. & Thimm, G. Systematic prediction of 70. Ceriotti, M., Tribello, G. A. & Parrinello, M. Simplifying the representation of crystal structures. Chem. Phys. Lett. 337, 36 (2001). complex free energy landscapes using sketch-map. Proc. Natl. Acad. Sci. USA 41. Winkler, B., Pickard, C. J., Milman, V., Klee, W. E. & Thimm, G. Prediction of 108, 13023 (2011). nanoporous sp -carbon framework structure by combining graph theory with quantum mechanics. Chem. Phys. Lett. 312, 536 (1999). Acknowledgements 42. Strong, R. T., Pickard, C. J., Milman, V., Thimm, G. & Winkler, B. Systematic E.A.E., R.J.N. and C.J.P. acknowledge ﬁnancial support from the Engineering and Phy- prediction of crystal structures: an application to sp -hybridized carbon sical Sciences Research Council of the UK [EP/J017639/1]. C.J.P. was also supported by a polymorphs. Phys. Rev. B 70, 045101 (2004). Royal Society Wolfson Research Merit Award. M.C. and A.A. acknowledge funding by 43. Baburin, I. A., Proserpio, D. M., Saleev, V. A. & Shipilova, A. V. From zeolite the European Research Council under the European Union’s Horizon 2020 research and nets to sp carbon allotropes: a topology-based multiscale theoretical study. innovation programme (grant agreement no. 677013-HBMAP). The calculations were Phys. Chem. Chem. Phys. 17, 1332 (2015). performed on the Cambridge High Performance Computing Service facility and the 44. Matsui, T., Hirata, M., Yagasaki, T., Matsumoto, M. & Tanaka, H. Archer facility of the UK’s national high-performance computing service (for which Hypothetical ultralow-density ice polymorphs. J. Chem. Phys. 147, 091101 access was obtained via the UKCP consortium [EP/P022596/1]). (2017). 45. Baerlocher, C. & McCusker, L. Database of zeolite structures. IZA Structure Commission http://www.iza-structure.org/databases/ (2017). Author contributions 46. Engel, E. A., Monserrat, B. & Needs, R. J. Anharmonic nuclear motion and the E.A.E., C.J.P. and R.J.N. led the structure search. E.A.E., A.A. and M.C. devised the relative stability of hexagonal and cubic ice. Phys. Rev. X 5, 021033 (2015). generalised convex-hull framework. A.A. and M.C. led the machine-learning analysis. All 47. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation authors contributed to the discussion and writing of the manuscript. made simple. Phys. Rev. Lett. 77, 3865 (1996). 48. Lee, K., Murray, E. D., Kong, L., Lundqvist, B. I. & Langreth, D. C. Higher- accuracy van der Waals density functional. Phys. Rev. B 82, 081101(R) (2010). Additional information 49. van Duin, A. C. T., Dasgupta, S., Lorant, F. & III, W. A. G. ReaxFF: a reactive Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- 018-04618-6. force ﬁeld for hydrocarbons. J. Phys. Chem. A 105, 9396 (2001). 50. del Rosso, L. et al. Reﬁned structure of metastable ice XVII from neutron diffraction measurements. J. Phys. Chem. C 120, 26955 (2016). Competing interests: The authors declare no competing interests. 51. von Stackelberg, M. & Müller, H. Zur struktur der gashydrate. Naturwissenschaften 38, 456 (1951). Reprints and permission information is available online at http://npg.nature.com/ 52. Claussen, W. Suggested structures of water in inert gas hydrates. J. Chem. reprintsandpermissions/ Phys. 19, 259 (1951). Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in 53. Pauling, L. & Marsh, R. The structure of chlorine hydrate. Proc. Natl. Acad. Sci. USA 38, 112 (1952). published maps and institutional afﬁliations. 54. Jeffrey, G. A. in Inclusion Compounds (eds Atwood, J. L., Davies, J. E. & MacNicol, D. D.) 135–185 (Academic Press, NY, 1984). 55. Rossi, M., Gasparotto, P. & Ceriotti, M. Anharmonic and quantum ﬂuctuations in molecular crystals: a ﬁrst-principles study of the stability of Open Access This article is licensed under a Creative Commons paracetamol. Phys. Rev. Lett. 117, 115702 (2016). Attribution 4.0 International License, which permits use, sharing, 56. Bartók, A. P., Kondor, R. & Csányi, G. On representing chemical adaptation, distribution and reproduction in any medium or format, as long as you give environments. Phys. Rev. B 87, 184115 (2013). appropriate credit to the original author(s) and the source, provide a link to the Creative 57. De, S., Bartók, A. P., Csányi, G. & Ceriotti, M. Comparing molecules and Commons license, and indicate if changes were made. The images or other third party solids across structural and alchemical space. Phys. Chem. Chem. Phys. 18, material in this article are included in the article’s Creative Commons license, unless 13754 (2016). indicated otherwise in a credit line to the material. If material is not included in the 58. Sun, W. et al. The thermodynamic scale of inorganic crystalline metastability. article’s Creative Commons license and your intended use is not permitted by statutory Sci. Adv. 2, e1600225 (2016). regulation or exceeds the permitted use, you will need to obtain permission directly from 59. Schölkopf, B., Smola, A. & Müller, K.-R. Nonlinear component analysis as a the copyright holder. To view a copy of this license, visit http://creativecommons.org/ kernel eigenvalue problem. Neural Comput. 10, 1299 (1998). licenses/by/4.0/. 60. Hermann, A., Ashcroft, N. W. & Hoffmann, R. High pressure ices. Proc. Natl. Acad. Sci. USA 109, 745 (2011). 61. Liu, Y. & Ojamäe, L. Clathrate ice sL: a new crystalline phase of ice with © The Author(s) 2018 ultralow density predicted by ﬁrst-principles phase diagram computations. Phys. Chem. Chem. Phys. 20, 8333 (2018). NATURE COMMUNICATIONS (2018) 9:2173 DOI: 10.1038/s41467-018-04618-6 www.nature.com/naturecommunications 7 | | |
Nature Communications – Springer Journals
Published: Jun 5, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.
Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.
All the latest content is available, no embargo periods.
“Hi guys, I cannot tell you how much I love this resource. Incredible. I really believe you've hit the nail on the head with this site in regards to solving the research-purchase issue.”Daniel C.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud
“I must say, @deepdyve is a fabulous solution to the independent researcher's problem of #access to #information.”@deepthiw
“My last article couldn't be possible without the platform @deepdyve that makes journal papers cheaper.”@JoseServera