Systematic comparison of crystalline and amorphous phases: Charting the landscape of water structures and transformations Fabio Pietrucci and Roman Martoňák Citation: The Journal of Chemical Physics 142, 104704 (2015); doi: 10.1063/1.4914138 View online: http://dx.doi.org/10.1063/1.4914138 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/10?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Rich structural phase diagram and thermoelectric properties of layered tellurides Mo1−x Nb x Te2 APL Mat. 3, 041514 (2015); 10.1063/1.4913967 Reflectance changes during shock-induced phase transformations in metals Rev. Sci. Instrum. 81, 065101 (2010); 10.1063/1.3430536 Plastic crystal phases of simple water models J. Chem. Phys. 130, 244504 (2009); 10.1063/1.3156856 Interaction potential for silicon carbide: A molecular dynamics study of elastic constants and vibrational density of states for crystalline and amorphous silicon carbide J. Appl. Phys. 101, 103515 (2007); 10.1063/1.2724570 Structural transformation of Sb x Se 100 − x thin films for phase change nonvolatile memory applications J. Appl. Phys. 98, 014904 (2005); 10.1063/1.1946197

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.157.188.52 On: Tue, 10 Mar 2015 13:39:26

THE JOURNAL OF CHEMICAL PHYSICS 142, 104704 (2015)

Systematic comparison of crystalline and amorphous phases: Charting the landscape of water structures and transformations Fabio Pietrucci1,a) and Roman Martoňák2 1 2

Sorbonne Universités, UPMC University Paris 6, UMR 7590, IMPMC, F-75005 Paris, France Department of Experimental Physics, Comenius University, Mlynská Dolina F2, 842 48 Bratislava, Slovakia

(Received 8 December 2014; accepted 24 February 2015; published online 10 March 2015) Systematically resolving different crystalline phases starting from the atomic positions, a mandatory step in algorithms for the prediction of structures or for the simulation of phase transitions, can be a non-trivial task. Extending to amorphous phases and liquids which lack the discrete symmetries, the problem becomes even more difficult, involving subtle topological differences at medium range that, however, are crucial to the physico-chemical and spectroscopic properties of the corresponding materials. Typically, system-tailored order parameters are devised, like global or local symmetry indicators, ring populations, etc. We show that a recently introduced metric provides a simple and general solution to this intricate problem. In particular, we demonstrate that a map can be traced displaying distances among water phases, including crystalline as well as amorphous states and the liquid, consistently with experimental knowledge in terms of phase diagram, structural features, and preparation routes. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4914138] I. INTRODUCTION

Discriminating different phases of matter starting from the knowledge of the atomic positions, by means of a suitably sensitive metric, is a fundamental problem that is relevant to different research fields. In particular, in the last years, there have been growing efforts in the areas of in silico crystal structure prediction and materials’ informatics, where techniques like genetic algorithms,1 particle-swarm optimization,2 and random searches3 are routinely employed to discover lowenergy atomic arrangements inside a specified crystal unit cell with specified number of atoms. An important ingredient of these approaches is a metric able to resolve different phases starting from the atomic positions, allowing to distinguish newly found structures from those already explored. This can be non-trivial and different solutions have been put forward, including fingerprints based on local topology and symmetry4,5 or based on pair correlation functions.6 A similar problem is faced in the simulation of phase transitions, e.g., employing techniques like umbrella sampling,7 metadynamics,8 or transition path sampling.9 In these cases, a suitable order parameter has to be identified that follows the progress of the structural transformation, allowing to trace a free energy landscape. Finding a general formulation of a metric that systematically distinguishes different phases of condensed matter can be non-trivial already in the case of crystalline materials, characterized by discrete symmetry and long-range order. The problem becomes even more difficult if one enlarges to consider also disordered systems such as amorphous materials and liquids that differ by their peculiar medium-range structural features. Taking, for instance, the case of structural transformations between different forms of water, a difficult system due to its numerous polymorphs a)Electronic mail: [email protected]

0021-9606/2015/142(10)/104704/7/$30.00

and anomalies, solutions tailored upon the chemistry of the substance and the specificities of the phases considered have been sought for, including Steinhardt’s local order parameters10–14 or populations of specific classes of rings.15 A conceptually different and more general approach has been successfully applied in metadynamics simulations based on cell parameters as order parameters;16,17 this method, however, is especially suited to structural phase transitions between crystals and only in presence of sizable changes in the shape or size of the unit cell. Recently, one of us introduced a new metric based on a very general formulation, obtained from the matrix of allto-all interatomic distances and invariant upon permutation of identical atoms or molecules, and showed that it allows tracking in an accurate way, the progress of activated processes in solution, including chemical reactions and ion pairing/dissociation.18 Here, we show that, based on this metric, it is possible to achieve the result of systematically assigning a structural distance between the crystalline and even amorphous phases of water, in an insightful way fully consistent with the available experimental knowledge. Furthermore, the metric allows tracking the detailed evolution of structural transitions among crystalline and amorphous water ices. Our approach is fully general, since we do not exploit any information about the system-specific symmetry or chemical bond topology, and in fact, we show at the end of this article that it is possible to pass in a seamless way from the study of water phases to the study of oxide glasses quenched from the melt. II. METHODS

We start from the definition of Permutation Invariant Vector (PIV) introduced in Ref. 18: such vector characterizes the structure of a configuration of N atoms in Cartesian space

142, 104704-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.157.188.52 On: Tue, 10 Mar 2015 13:39:26

104704-2

F. Pietrucci and R. Martoňák

{Ri }i=1, .., N and is constructed as     v = sort C(∆Ri j ) = sort C(|Ri − R j |) , i > j,

J. Chem. Phys. 142, 104704 (2015)

(1)

where C is a switching function decaying from one to  −1 zero (here, we adopt C(∆Ri j ) = 1 + e(∆R i j −ρ)/λ ) aimed at focusing on a specified range of interatomic distances. In the present work, we tested functions C with ranges from 1–3.5 Å up to 1–6 Å. The components of the N(N − 1)/2 dimensional vector v are sorted in ascending order within each set corresponding to the same pair of elements (e.g., O–O, O–H, and H–H for water). The latter operation ensures invariance of the vector upon permutation of identical atoms: this is a crucial property in view of an efficient representation and exploration of the configuration space, especially when dealing with disordered non-crystalline systems (see also Ref. 19 for a related approach focused on isolated molecules and atomic clusters). Two atomic configurations A and B can be then compared by means of the Euclidean distance between the corresponding PIVs v A and v B,     1 N (N −1)
2 1 2  vkA − vkB , (2) D AB = |v A − v B | = N k=1 where, at variance with Ref. 18, we introduced the system size scaling factor 1/N to facilitate the comparison of absolute values of D among systems with different numbers of atoms. To perform all calculations presented here, we developed the freely downloadable software tool piv_clustering version 1.3.20 Starting from the PIV distances between all pairs of structures in a given set, for the sake of visualization and clarity, a two-dimensional map can be constructed employing the following Monte Carlo optimization algorithm. A representative point in a plane is assigned to each structure, initially with random position (0 < x < 1, 0 < y < 1). Next, a cost function is defined as  U= (Di j − d i j )2, (3) i< j

where the sum runs over all pairs of structures and each harmonic term grows with the difference between the PIV distance Di j and the in-plane distance d i j . In the case of liquid water and amorphous ices in Fig. 3, a set of 50 configurations is considered for each disordered phase but for the sake of clarity, a single point is shown for each phase, with Di j defined as the average distance between the configurations of phase i and those of phase j. U is minimized through a sequence of moves, each move consisting in the random displacement (maximum amplitude 0.0005) of one point in the plane. Moves are accepted or rejected according to a Metropolis criterion at thermal energy kT = 1 × 10−7 (in U units) until obtaining a Pearson’s correlation coefficient = 0.999 between PIV distances and in-plane distances (in the case of water phases in Fig. 3, 105 moves were sufficient). We verified that the resulting map is robust with respect to variations in the choice of the Monte Carlo parameters (initial positions, temperature, and move amplitude): when the correlation coefficient reaches a value of 0.999, the very same

map is repeatedly obtained up to irrelevant rotations. The map conveys in a direct way the information about the topological relationships in the set of structures, in the same spirit as scatter plots of crystal structures in Ref. 21, what could also be obtained employing more complex dimensionality reduction algorithms like sketch maps.22 We remark that the two axes of the plane are arbitrary and, in principle, they do not bear a physical meaning.

III. RESULTS A. Analysis of water phases

Water has a very complex phase diagram: we consider here the liquid together with several molecular crystalline polymorphs of ice (Ic , Ih , II, III, VI, VII, XII)23 as well as various amorphous phases including the three wellknown ones—Low Density Amorphous (LDA), High Density Amorphous (HDA), and Very High Density Amorphous (VHDA) ices.24 All these phases of water received extensive experimental and theoretical attention and are characterized by different symmetries and topologies of the hydrogen bond network. In order to compute the distances D among all phases, we generated a sequence of atomic configurations, each one including 360 water molecules in a periodically repeated supercell. In the case of ices II, III, VI, VII, and XII, one configuration for each phase is included in the sequence, which is taken from the crystallographic databases (adding a random occupancy of proton sites in ice VII and considering only the oxygen positions in ices III, VI, and XII). In the case of ices Ic and Ih , 20 configurations for each phase are included, differing by their random population of hydrogen sites. Fifty configurations are included for each of the amorphous ices, sampled from 0.5 ns-long molecular dynamics (MD) simulations reported in Ref. 25, in the isobaric-isothermal N pT ensemble at p = 0 and T = 80 K. Water was modeled by employing the TIP4P interatomic potential26 that is able to reproduce quite accurately the topology of the phase diagram of water.27 Amorphous ices were prepared by the protocol described in Ref. 25 and decompressed to zero pressure and are labeled by the pressure at which the thermal annealing was performed. The HDA ice represents the as-prepared outcome of pressure-induced amorphization of Ih ice at T = 80 K decompressed from p = 15 kbar to p = 0 without any thermal annealing (in Ref. 25, this ice was denoted as HDA’). Amorphous ice annealed at p = 19.5 kbar and decompressed to p = 0 is labeled as VHDA. Finally, in the case of liquid water, a new 50 ns-long N pT MD simulation at ambient conditions was performed with the same settings as the previous ones, sampling 50 configurations from the last 40 ns of trajectory. Next, we analyzed the sequence of all configurations described above by computing the matrix of all-to-all distances D, employing a function C that focuses between ≈1 and ≈4.5 Å (with ρ = 2.6 Å, λ = 0.6 Å). Since including all atoms (oxygen and hydrogen) in the definition of PIV, we found very similar results as including only oxygen; we focus our discussion on the latter case (for a full comparison see the supplementary material28). The first important result

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.157.188.52 On: Tue, 10 Mar 2015 13:39:26

104704-3

F. Pietrucci and R. Martoňák

J. Chem. Phys. 142, 104704 (2015)

considered. Due to the limited range of the switching function C, the components become very small beyond about the 20th neighbor of each atom. Clearly ices Ih and Ic are highly symmetrical and very similar to each other, with a first shell of four neighbors and a second shell of 12 atoms, their PIVs differing mainly for an extra atom close to the second shell that is present only in Ih . Ices II, III, and XII feature a more graded set of distances in the first as well as in the second shell. Ice VI displays a second shell with sizeably shorter interatomic distances with respect to the previous crystals, suggesting a higher density, whereas ice VII further continues in this trend and is peculiar in having a first shell composed by eight atoms instead of four. The three amorphous forms LDA,

FIG. 1. Distribution of PIV distances between pairs of water phases (only oxygen atoms are included in the definition).

is that, despite of the thermal fluctuations, in the case of amorphous states and liquid water, the PIV-based metric is able to discriminate whether two configurations belong to the same phase or to different ones. As shown in Fig. 1, the distances between configurations belonging to the same phase are systematically smaller (<0.02) than the distances between configurations belonging to different phases (>0.09), with a clear gap in between. The discriminating power of the PIV metric is confirmed by a further benchmark: we generated four independent new samples of HDA ice by repeating the procedure in Ref. 25, each time starting from Ih ice with a different set of randomized initial velocities at T = 80 K, compressing up to 15 kbar within 20 ns and decompressing to zero pressure within 55 ns, thus ending up in four different geometric realizations of the amorphous material. For each new sample, 50 configurations were extracted from a 0.5 ns-long zero-pressure simulation, as in the case of the original sample from Ref. 25. The PIV distances between configurations belonging to different HDA samples (average = 0.0087, standard deviation = 0.0023) turn out to be slightly larger than those between configurations belonging to the same sample (average = 0.0070, standard deviation = 0.0015), but still well below 0.02, thus out of the gap region in Fig. 1. We remark that the resolving power of the PIV metric is not critically reduced neither by time correlations within the atomic trajectory of a given phase nor by mere density variations, as discussed in detail in the supplementary material.28 Fig. 2 displays the PIV components vi , in decreasing order, for one representative structure of each water phase

FIG. 2. Graphical representation of the PIV components v i associated to the different water phases considered. On the horizontal axis, the component index i is divided by the number of atoms N to facilitate the identification of shells of neighbors.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.157.188.52 On: Tue, 10 Mar 2015 13:39:26

104704-4

F. Pietrucci and R. Martoňák

FIG. 3. (a) Optimized two-dimensional map (employing arbitrary axes) that accurately reproduces the PIV distances between water phases: the size of spots corresponding to liquid and LDA, HDA and VHDA phases matches the spread of PIV distances within the ensemble of configurations (Fig. 1); (b) magnification of the region of amorphous ices. The red points are labeled by the pressure (in kbar) at which the corresponding amorphous structures were obtained. (c) Experimental phase diagram of water.29,30

HDA, and VHDA appear at first sight as smoothed versions of the crystalline polymorphs, with a markedly decreasing distance of atoms beyond the first shell for increasing density. Finally, the inflection point that in amorphous ices suggests the transition from the first to the second shell is barely visible in liquid water. For the sake of comparison, we plotted also the PIV profile of a random configuration of hard spheres of radius 1.25 Å (half the minimal O–O distance in water) at liquid water density (volume fraction = 0.27) and found that it differs from that of liquid water mostly in the first shell. A central result of this work is represented by Fig. 3, where we charted the structural landscape of the different water phases. This two-dimensional “topological map” despite having arbitrary axes represents in a quantitatively faithful manner

J. Chem. Phys. 142, 104704 (2015)

the true PIV distances between phases (Pearson’s correlation coefficient = 0.99928) and is constructed through the simple Monte Carlo algorithm described in Sec. II. It is remarkable that the PIV distances are so well representable in 2D. Comparison between this topological map and the canonical (p,T) phase diagram of water29,30 shows a number of interesting similarities (note that our data refer to ambient pressure for all phases except for ice II obtained at 3 kbar,31 ice III at 2.8 kbar,32 ice VII at 24 kbar,33 and ice XII at 5 kbar29). To begin with, the ranking of topological distances III < II < XII < VI < VII with respect to liquid water (at normal conditions) qualitatively reflects the distance from ambient-conditions liquid in the phase diagram along the pressure axis. Ices Ic (metastable) and Ih , close to each other in the map, are known to be structurally very similar to each other and also adjacent in the phase diagram. Remarkably, in the topological map, amorphous ices are located between liquid water and crystalline ices, correctly suggesting that they could be reached “on the way” to freezing toward the crystals, what is achieved in experiments at very high cooling rates.34 This localization along the “melting line” is also consistent with the experimental observation that upon heating at different pressures, amorphous ices display a recrystallization behavior that follows closely the crystallization behavior of liquid water.35 More in detail, LDA appears on the way between liquid and ice I, whereas HDA and VHDA are closer to ice II. This is consistent with experimental preparation routes: upon heating, LDA crystallizes into ice Ic ,34,36 and the two forms were also found to have similar features in vibrational spectrum.37 HDA ice together with ices II and VI (approximately equidistant from ice VII in the topological map) is known to be able to convert to ice VII upon compression to 4 GPa at low temperature,38,39 whereas ice Ih (much more distant from ice VII in the topological map) does not show this transition. Also VHDA ice has been reported to convert into VII at 3-4 GPa,35 again in agreement with the proximity in the topological map. Experiments and simulations suggest that VHDA has some structural similarity with ice VII since both contain interstitial non-bonded molecules close to the first shell, but disordered VHDA ice lacks the interpenetrated hydrogen networks present in ice VII.40–42 The proximity of HDA ice to ice II also suggests that the pressureinduced amorphization of ice Ih resulting instead of ice II in HDA can be regarded as an attempt to create ice II which is, however, frustrated due to kinetic reasons related to low atomic mobility at low temperatures. On the other hand, we are not aware of experiments directly connecting ice III with LDA or HDA, despite their proximity in the map. In this case, it is possible that the very limited stability range of ice III in the phase diagram makes it arduous to realize a direct transformation without passing through other crystalline phases of broader stability (Fig. 3(c)). In summary, the map not only displays distances between phases that are consistent with available experimental and theoretical structural information, but also in several cases, proximity in the map corresponds to the existence of a viable preparation route, even if there are also cases where a direct connection has not (yet) been reported. As a further evidence, note how the structures of amorphous ices obtained at increasing pressures (from 0-LDA to 19.5 kbar-VHDA, each one

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.157.188.52 On: Tue, 10 Mar 2015 13:39:26

104704-5

F. Pietrucci and R. Martoňák

J. Chem. Phys. 142, 104704 (2015)

then relaxed to ambient pressure) describe a path in the map that leads from LDA through HDA until VHDA polyamorphs (Fig. 3(b)). Here, we also observe that the distances between the amorphous phases formed at equidistant pressure intervals are first small at pressures below 1.5 kbar, then grow very fast between 1.5 and 3 kbar as the system transforms from LDA to HDA megabasin and upon further compression, decrease again as the limiting VHDA form is approached. This suggests that the PIV metric correctly captures the salient features of the structural evolution of a complex material like water, as further confirmed by the analysis of phase transitions that follows. It is also tempting to approximately interpret the two dimensions of the topological map (after suitable rotation, they can always be oriented as in Fig. 3(a)) as density and order, the former increasing in the horizontal direction and the latter decreasing in the vertical direction. Liquid water then represents the most disordered phase in this diagram which is plausible. To our knowledge, it is the first time that a comprehensive and insightful map like the one presented here can be drawn for a number of water phases starting from a single metric based on plain Cartesian coordinates. B. Structural transitions in crystalline and amorphous water

We explicitly tested the applicability of the PIV approach also to reactive trajectories including crystalline as well as amorphous phases. It is well-known that upon compression to 12 kbar at liquid nitrogen temperature, the pressure-induced amorphization from hexagonal ice Ih to HDA occurs in abrupt manner similar to a first-order phase transition.43,44 The same is true for the compression of LDA to 6 kbar.44 Both transformations were studied by constant-pressure MD in Refs. 25 and 45, where they were also found to be discontinuous in agreement with experiment. We applied the PIV approach first to the amorphization of Ih ice which was compressed from 0 to 13.5 kbar in steps of 1.5 kbar and at each pressure, a 5.5 ns MD run was performed (for further details, see Ref. 25). At 13.5 kbar, the Ih structure collapsed and transformed in the HDA amorphous ice. We analyzed the first 1 ns of this trajectory and the results are shown in Fig. 4. In panel (a), we see that the PIV metric clearly identifies five structurally homogeneous parts of the trajectory. Comparing these parts with the time evolution of density (b) and total energy (c), we see that they correlate with the plateaus of the latter quantities especially at the later stages of the transformation. For comparison, we show in Fig. 4(d) also the radial distribution functions averaged over the structurally homogeneous regions. This clearly confirms that the regions distinguished by the PIV metric indeed exhibit different short and medium range orders. In Fig. 5, we show the application of the metric to the LDA → HDA transformation which represents a paradigmatic pressure-induced transformation between two amorphous phases. The simulation protocol was the same as above. The transition takes place in two stages, starting at 7.5 kbar and completing at 9 kbar. In Fig. 5, we analyze the whole compression up to 22.5 kbar. Here, we see that even in case of transformation between amorphous phases, the PIV metric

FIG. 4. Pressure induced amorphization Ih → HDA at T = 80 K and p = 13.5 kbar. (a) PIV distance between configurations; time evolution (b) of the density and (c) of the total energy during the simulation; (d) O–O pair distribution function averaged over the stationary regions (metastable states) of panel (a).

is a much more convenient indicator of structural change than energy and density. The latter quantities grow upon each pressure increase, even when the compression is purely elastic and not accompanied by any change of the network topology. On the other hand, the PIV metric clearly identifies the points at 7.5 and 9 kbar where the inelastic densification takes place due to network reconstruction. In Fig. 5(a), one can see that the additional densification of the structure upon

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.157.188.52 On: Tue, 10 Mar 2015 13:39:26

104704-6

F. Pietrucci and R. Martoňák

J. Chem. Phys. 142, 104704 (2015)

form of a glass (Tg ≈ 580 K, Tm ≈ 1000 K), characterized by a covalent bond network of complicated topology, with tellurium atoms adopting coordinations between 3 and 5 and oxygen atoms between 1 and 2. We analyzed the ab initio MD trajectories in Ref. 46, where a 96-atoms liquid model was quenched from 2400 K down to 300 K obtaining a glass model that reproduces the neutron diffraction47,48 and infrared/Raman spectra49 of the real glass. This system represents a difficult test case for the PIV metric due to the similarities between the glass and the melt in terms of coordination shells46 and due to the volume being kept artificially constant throughout, thus suppressing density variations. The PIV metric combined with a structural cluster analysis (performed as in Ref. 18) is able to automatically group together structures generated at the same temperature, separating with high accuracy the glass from the melt at 1000 K and at 2400 K. As shown in Fig. 6, the typical PIV distance (computed including tellurium and oxygen atoms) between an atomic configuration of the glass at 300 K and a configuration of the melt at 1000 K is larger than the typical distance between two configurations at the same temperature; however, at variance with the case of water, the distributions overlap sizably, due to the aforementioned structural similarity between the two phases (the liquid at 2400 K appears instead very different from the other two). To test the resolving power of the PIV metric, we assembled a sequence containing 300 atomic configurations from microcanonical NVT ab initio MD simulations, 100 for each of the temperatures considered (300, 1000, and 2400 K), and following Ref. 18, we performed a cluster analysis employing the robust k-medoids algorithm:50 1. A specified number n of cluster centers is initially selected among all configurations based on a stochastic procedure (k-means++ initialization51) and clusters are formed by Voronoi partitioning, i.e., by assigning each configuration to its closest center. 2. The n centers are updated in order to minimize the sum of distances between the center and the members of each cluster, and Voronoi partitioning is again applied.

FIG. 5. Pressure induced amorphization LDA → HDA at T = 80 K: a series of constant pressure MD trajectories are analyzed, between 1.5 and 22.5 kbar. (a) PIV distance between configurations; time evolution (b) of the density and (c) of the total energy during the simulation; (d) O–O pair distribution function averaged at different pressures.

Update step (2) is repeated until the centers do not change anymore. To improve statistics, we repeated 10 times the cluster analysis, and to assess the robustness of the PIV metric,

further compression to 22.5 kbar is accompanied only by minor changes in the metric. This reveals that no more major structure changes are going on because of kinetic barriers at low temperature (see also the radial distribution function in Fig. 5(d)). C. TeO2 glass and liquid

Finally, we turn to the question of how general is the applicability of our new approach by analysing a different system, tellurium oxide (TeO2). This material can assume the

FIG. 6. Distribution of PIV distances between pairs of TeO2 structures at the same temperature (e.g., 300 K vs 300 K), compared with the distribution between pairs at different temperatures (e.g., 300 K vs 1000 K). A C function with ρ = 2.4 Å, λ = 0.4 Å is employed.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.157.188.52 On: Tue, 10 Mar 2015 13:39:26

104704-7

F. Pietrucci and R. Martoňák

we considered, as for water, different C functions with ranges from 1–3.5 Å up to 1–6 Å. When the configuration space is partitioned into two clusters, the first one contains a mix of structures at 300 and 1000 K and the second one only structures at 2400 K, as it could be expected from Fig. 6. Partitioning into three clusters, on an average, each of them has 93% of its members belonging to the same temperature (either glass at 300 K, or liquid at 1000 K, or liquid at 2400 K), and partitioning into four clusters, the figure rises to 98%. This indicates that the PIV metric is still capable of automatically discriminating with high accuracy the different sets of structures, with a weak dependence of the results (variations in populations of about 2%) upon the range of the function C.

IV. CONCLUSIONS

Employing without any fine tuning the metric based on PIVs, it is possible to resolve well all the distinct (ordered and disordered) phases we considered, in the case of water as well as in the case of tellurium dioxide, and it is also possible to analyze in a detailed and insightful way structural transformations like the amorphization of water ice or the polyamorphic transition between LDA and HDA. The computational approach we propose has the advantage of being simple and general, and it can be exploited in several ways. In combination with algorithms for structure prediction1–3 and with databases for materials’ informatics, it can provide an important additional insight in the form of a map suggesting the viable preparation routes connecting the structures explored. In combination with simulation techniques for the study of phase transitions, it can form the basis of new order parameters (already under study) able to address the dynamic mechanisms connecting together crystalline as well as liquid and amorphous structures, hopefully also leading to a better understanding of glasses and the glass transition. ACKNOWLEDGMENTS

We gratefully acknowledge Wanda Andreoni, Marco Saitta, Stefan Klotz, and Marián Fecko for insightful discussions and suggestions. R.M. was supported by the Slovak Research and Development Agency under Contract Nos. APVV-0558-10 and APVV-0108-11. R.M. also acknowledges support from CECAM during a collaborative visit. 1C.

Glass, A. Oganov, and N. Hansen, Comput. Phys. Commun. 175, 713 (2006). 2Y. Wang, J. Lv, L. Zhu, and Y. Ma, Phys. Rev. B 82, 094116 (2010). 3C. Pickard and R. Needs, J. Phys.: Condens. Matter 23, 053201 (2011). 4V. A. Blatov, IUCr CompComm Newsletter 7, 4 (2006). 5R. Hundt, J. Schon, and M. Jansen, J. Appl. Crystallogr. 39, 6 (2006). 6A. Oganov and M. Valle, J. Chem. Phys. 130, 104504 (2009). 7G. M. Torrie and J. P. Valleau, J. Comput. Phys. 23, 187 (1977). 8A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A. 99, 12562 (2002). 9P. Bolhuis, D. Chandler, C. Dellago, and P. Geissler, Annu. Rev. Phys. Chem. 53, 291 (2002).

J. Chem. Phys. 142, 104704 (2015) 10P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983). 11R.

Radhakrishnan and B. L. Trout, J. Am. Chem. Soc. 125, 7743 (2003). Quigley and P. M. Rodger, J. Chem. Phys. 128, 154518 (2008). 13J. C. Palmer, R. Car, and P. G. Debenedetti, Faraday Discuss. 167, 77 (2013). 14D. T. Limmer and D. Chandler, Proc. Natl. Acad. Sci. U. S. A. 111, 9413 (2014). 15D. Donadio, P. Raiteri, and M. Parrinello, J. Phys. Chem. B 109, 5421 (2005). 16R. Martoˇ nák, A. Laio, and M. Parrinello, Phys. Rev. Lett. 90, 075503 (2003). 17R. Martoˇ nák, D. Donadio, A. R. Oganov, and M. Parrinello, Nat. Mater. 5, 623 (2006). 18G. A. Gallet and F. Pietrucci, J. Chem. Phys. 139, 074101 (2013). 19F. Pietrucci and W. Andreoni, Phys. Rev. Lett. 107, 085504 (2011). 20See http://www.sourceforge.net/projects/pivclustering. 21M. Valle and A. R. Oganov, Acta Crystallogr., Sect. A: Found. Crystallogr. 66, 507 (2010). 22M. Ceriotti, G. A. Tribello, and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A. 108, 13023 (2011). 23We consider here all molecular ices except for the metastable ice IV and stable ice V for which the number of molecules in the unit cell (16 and 28, respectively) is not commensurate with the amorphous samples, thus precluding the direct comparison. 24T. Loerting and N. Giovambattista, J. Phys.: Condens. Matter 18, R919 (2006). 25R. Martoˇ nák, D. Donadio, and M. Parrinello, J. Chem. Phys. 122, 134501 (2005). 26W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). 27E. Sanz, C. Vega, J. Abascal, and L. MacDowell, Phys. Rev. Lett. 92, 255701 (2004). 28See supplementary material at http://dx.doi.org/10.1063/1.4914138 for computational details and further analyses. 29C. Lobban, J. Finney, and W. Kuhs, Nature 391, 268 (1998). 30S. Fanetti, M. Pagliai, M. Citroni, A. Lapini, S. Scandolo, R. Righini, and R. Bini, J. Phys. Chem. Lett. 5, 3804 (2014). 31B. Kamb, Acta Crystallogr. 17, 1437 (1964). 32J. Londono, W. Kuhs, and J. Finney, J. Chem. Phys. 98, 4878 (1993). 33W. Kuhs, J. Finney, C. Vettier, and D. Bliss, J. Chem. Phys. 81, 3612 (1984). 34G. Johari, A. Hallbrucker, and E. Mayer, Nature 330, 552 (1987). 35S. Klotz, G. Hamel, J. S. Loveday, R. J. Nelmes, and M. Guthrie, Z. Kristallogr. - Cryst. Mater. 218, 117 (2003). 36A. Hallbrucker and E. Mayer, J. Phys. Chem. 91, 503 (1987). 37H. Schober, M. M. Koza, A. Tölle, C. Masciovecchio, F. Sette, and F. Fujara, Phys. Rev. Lett. 85, 4100 (2000). 38J. L. Finney, C. Lobban, S. Klotz, and J. M. Besson, in Collected Abstracts of the XVII Congress and General Assembly, edited by J. F. Griffin (IUCr, Seattle, 1996), Vol. C-534; J. L. Finney, C. Lobban, S. Klotz, J. M. Besson, and G. Hamel, in ISIS 96—Rutherford Appleton Laboratory Report RAL 96050, edited by A. D. Taylor (Chilton, UK, 1996), Vol. A17. 39S. Klotz, J. Besson, G. Hamel, R. Nelmes, J. Loveday, and W. Marshall, Nature 398, 681 (1999). 40A. M. Saitta, T. Strässle, G. Rousse, G. Hamel, S. Klotz, R. J. Nelmes, and J. S. Loveday, J. Chem. Phys. 121, 8430 (2004). 41A. M. Saitta, T. Strässle, and S. Klotz, J. Phys. Chem. B 110, 3595 (2006). 42A. Saitta, T. Strässle, and S. Klotz, Europhys. Lett. 74, 445 (2006). 43O. Mishima, L. D. Calvert, and E. Whalley, Nature (London) 310, 393 (1984). 44O. Mishima, L. D. Calvert, and E. Whalley, Nature (London) 314, 76 (1985). 45R. Martoˇ nák, D. Donadio, and M. Parrinello, Phys. Rev. Lett. 92, 225702 (2004). 46F. Pietrucci, S. Caravati, and M. Bernasconi, Phys. Rev. B 78, 064203 (2008). 47H. Niida, T. Uchino, J. Jin, S.-H. Kim, T. Fukunaga, and T. Yoko, J. Chem. Phys. 114, 459 (2001). 48E. R. Barney, A. C. Hannon, D. Holland, N. Umesaki, M. Tatsumisago, R. G. Orman, and S. Feller, J. Phys. Chem. Lett. 4, 2312 (2013). 49O. Noguera, T. Merle-Méjean, A. Mirgorodsky, M. Smirnov, P. Thomas, and J.-C. Champarnaud-Mesjard, J. Non-Cryst. Solids 330, 50 (2003). 50H.-S. Park and C.-H. Jun, Expert Syst. Appl. 36, 3336 (2009). 51D. Arthur and S. Vassilvitskii, Proceedings of the Eighteenth Annual ACM–SIAM Symposium on Discrete Algorithms, SODA ’07 (Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2007), pp. 1027–1035. 12D.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.157.188.52 On: Tue, 10 Mar 2015 13:39:26