- Email: [email protected]
Density and porosity of amorphous water ice by DFT methods
Chemical Physics Letters 745 (2020) 137222
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Research paper
Density and porosity of amorphous water ice by DFT methods a,⁎
b
Pedro C. Gómez , Miguel Á. Satorre Aznar , Rafael Escribano a b c
T
c
Departamento de Química Física, Facultad de CC. Químicas, Universidad Complutense, and Unidad Asociada Physical Chemistry UCM-CSIC, 28040 Madrid, Spain Departamento de Física Aplicada, EPS Alcoi, Centro de Tecnologías Físicas, Universitat Politècnica de València, Placeta Ferrándiz-Carbonell s/n, 03801 Alcoi, Spain Instituto de Estructura de la Materia, IEM-CSIC, and Unidad Asociada Physical Chemistry UCM-CSIC, Serrano 123, 28006 Madrid, Spain
H I GH L IG H T S
water ice studied at Density Functional Theory level. • Amorphous (LDA) and high-density (HDA) water ices studied at 40 K and 140 K. • Low-density with experimental results on density, porosity and internal structure. • Agreement • Possibility to apply this method to study i.e. astrophysical ices of other species.
A R T I C LE I N FO
A B S T R A C T
Keywords: Amorphous ices Density Porosity DFT methods
A theoretical method is developed to study amorphous ices, based on solid state calculations at DFT level. The method is applied to amorphous water ice of a large range of densities, between 0.6 g cm−3 and 1.6 g cm−3 to represent low-density and high-density amorphous ices, at 40 K and 140 K. The most stable forms appear for densities in good agreement with experimental data, ρ ~ 0.9–1.0 g cm−3 and ρ ~ 1.1–1.2 g cm−3. The method allows evaluating the porosity through the intrinsic density, providing a possible way to estimate the porosity of astrophysical ices.
1. Introduction
2. Methodology
There may be already dozens of high value scientific papers dealing with amorphous water ice (usually denominated ASW) and its physical properties. These include, to quote just a few, density [1,2], which comprises low density, high density and very high density amorphous ice (LDA, HDA and VHDA, respectively); internal structure including Hbonding [3–8]; energetics and stability [4,9–11]; experimental production [9,12,13]; or theoretical models [14–18]. Within the latter topic, the only justification for one more paper may arise if that paper conveys some new methodology allowing the estimation of not sufficiently understood properties of the ice. Our aim in this investigation is thus to present a high-level theoretical methodology at moderate computational cost, test its validity on some well-known examples, and apply it to study ice porosity, which is clearly an underdeveloped subject, especially for astrophysical ices. We can also speculate on the application of this method to deal with an extremely interesting issue, that of the formation of a crystal.
The approach used in this investigation is based in the consideration of an amorphous solid as a 3-dimensional repetitive unit cell where a number of molecules are located in an unsystematic way. The disorderly arrangement is achieved by unlocking the symmetry of a crystalline structure through a fast force-field molecular dynamics using the Amorphous Cell modulus of the Materials Studio (MS) package. In our models, we use a cubic unit cell containing 12 water molecules. The density of the sample is thus determined by the length of the cell side a. This amorphization process has been carried out at two temperatures, Tamorph = 40 K and 140 K, and at a broad range of density values, for ρ between 0.6 g cm−3 and 1.6 g cm−3 (cubic cell side 8.426 Å and 6.073 Å, respectively); this range of densities spans a wider interval than normally achieved experimentally. Bearing in mind that the most abundant OeO distances in amorphous samples have been measured at ~2.7 and 4.3 Å [19,20], the size of our samples guarantees in every case that the most important molecular interactions are enclosed. The size of the cell and the number of molecules within, are very important factors. We have chosen 12 molecules inside the cell, which
⁎
Corresponding author. E-mail address: [email protected] (P.C. Gómez).
https://doi.org/10.1016/j.cplett.2020.137222 Received 21 December 2019; Received in revised form 10 February 2020; Accepted 12 February 2020 Available online 15 February 2020 0009-2614/ © 2020 Elsevier B.V. All rights reserved.
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
is not a large number but still allows describing intermolecular interactions, like H-bonding and van der Waals forces, using a fairly high level of quantum chemistry theory. Samples with hundreds of water molecules have been studied by much less demanding methods, like force-field dynamics models [21–23]. Quantum chemistry models of higher theoretical level are sometimes used for larger number of molecules, 32, 64 or even higher [24,25], but they demand the use of large computation facilities. With our model, we detect the formation of strong H-bonds among two or three molecules within a cell, which often continue in a chain-like form with molecules in the adjacent cells for higher density cases. In our samples, we have observed water molecules with usual coordination numbers from two to four, but including coordination number five as well in some cases [7]. These facts support the choice of 12 molecules per cell to be large enough to describe the main interactions within the ice. Moreover, we can apply our method to big amount of samples to cover an ample range of densities. The total number of samples that we have considered amounts to 230 at fixed cell geometry, and the same number for free cell calculations, discussed below. Our models do not aim to provide an accurate replica of the experiments, as indeed no model can be, but allow studying interesting physical properties. The validity of this methodology to predict new information must be tested by comparing to available experimental results, as we are doing in this investigation. Geometrical structures are optimized looking for a minimum in their potential energy surface (PES). This process is done at Density Functional Theory (DFT) level by means of the CASTEP [26] modulus of the MS package. CASTEP uses pseudopotentials and plane-wave basis sets. The k-point is modified along the increasing accuracy of our calculations. It ranges from (1 × 1 × 1), with k-point separation of 0.1 Å−1 at the initial stages to (2 × 2 × 2), with separation of 0.07 Å−1 at the final optimization process. We have chosen the Generalized Gradient Approximation (GGA) functional by Perdew-Burke-Ernzerhof (PBE) including Grimme’s D2 dispersion correction, which we have found to be very important when dealing with water molecules [16], and a plane-wave basis set cutoff of 830 eV.
Fig. 1. Electronic energy vs density for the set of selected densities. Each energy corresponds to the average of 10 optimized structures generated at 40 K (squares) and 140 K (circles). Zero energy for the drawing is set at −5564.4 eV, with 140 K values offset by −1.0 eV. The dotted line is only drawn to guide the eye.
The graph for 140 K has two local low values, at ~1.0 and ~1.3 g cm−3. There is some likelihood between both graphs, but the distribution of energies around the minima is different, and the high density local minimum for T = 40 K at ρ ~ 1.5 cm−3, is missing at T = 140 K. In the latter case, the structures in the 0.9–1.1 g cm−3 range are energetically close, whereas the local low at 1.2–1.3 g cm−3 is more marked. The interest of these results arises from the assumed inverse relationship between electronic energy and stability of the sample. Thus, we can expect the most stable amorphous structures to be formed at densities close to the low values shown in Fig. 1. Amorphous experimental ices are formed by isles or regions of water aggregates of different densities which give an overall total density of the sample. We can assume that the densities that look more favorable in our calculations will be dominant within the island structure of the experimental ice. The experimental measurements give densities of 1.1–1.2 g cm−3 for HDA grown at T < 40 K [3,28,29], which can be as high as 1.31 g cm−3 at ice-pressurized experiments.[9] On the other hand, for higher temperature experiments, between 70 K and up to near 130 K, LDA is generated with density close to that of cubic ice, 0.94 g cm−3 [20,30]. We can see that our calculations predict the most favorable structures to be formed with densities near the experimental observations, within the expected uncertainties.
3. Results 3.1. Average energy vs density For each Tamorph and ρ we generate 10 amorphous samples which are relaxed to optimize their structure. We choose 10 samples to ensure diversity of initial structures, and after the optimization, we average their energy and use this result, Eaverage, as a representation of the Tamorph vs ρ set. Fig. 1 collects the corresponding results. A full table with the energy and density of each relaxed structure is given in Table S1 of the Supplementary Material. The temperature of Tamorph = 40 K was chosen because it is considered to be the maximum temperature for the formation of HDA water ice.[27] In turn, Tamorph = 140 K was selected because it is close to the formation of crystalline cubic water ice Ic. These two temperatures are sufficiently far apart to permit a better appreciation of differences in the results.
3.2. Cell release It is also interesting to check the effects of releasing the size and shape of the cell together with the molecular arrangement. In this process, the cell will shrink or expand, depending on the internal stress of the structure, searching for a lower energy minimum, and yielding a new volume and therefore a new density. Thus, we have obtained another set of data, which now includes free-cell energies, Efree-cell, and released densities, ρfree-cell. The results are shown in Fig. 2 for Tamorph = 40 K (top) and 140 K (bottom). As above, average values over 10 structures are plotted in every case. Arrows join (Eaverage, ρ) of fixed-cell optimized structures, as in the previous section (squares), with (Efree-cell, ρfree-cell) of free-cell reoptimized structures (circles). The energy always decreases, sometimes as much as 1 eV, and the pressure, understood as the internal stress created in the cell during the
Tamorph = 40 K Fig. 1 indicates that there exist some values of the density which are energetically more favorable than their neighbours. At 40 K, they are located at ρ ~ 1.1 g cm−3, ~1.3 g cm−3, and ρ ~ 1.5 cm−3. The last one is clearly less promising than those in the lower density region. The Eaverage values in the 1.2–1.4 g cm−3 range are fairly close (see also Table S1 for the precise values), indicating that structures with these density values are almost isoenergetic, with the one at ρ ~ 1.1 g cm−3 being only 0.06 eV higher (1.4 kcal mol−1).
Tamorph = 140 K 2
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
Finally, in panel 4, we present a composition of our PCF for ρ = 1.3 g cm−3 and the trace obtained from the values tabulated by Narten et al. [20] To facilitate the comparison, Table 1 collects the measured values of OeO distances from each trace, and their assignment in terms of OeO characteristics. The most prominent peak, i.e. the most frequently occurring OeO distance, always corresponds to the separation between two O atoms linked via H-bond, O⋯HeO. The peak at ~4.5 Å, interpreted as OeO distance between non-bound OeO pairs, is stronger in less-crowded LDA samples, and weakens at higher densities, where water molecules are closer to each other, developing a peak at ~3.5 Å, ascribed to an interstitial O atom linked to other molecules by weak dispersion forces, giving a 5-coordination number O atom [14,19,20]. Finally, the last outstanding peak, at ~6.7 Å, corresponds to distant OeO pairs, which in our calculations are found between atoms in adjacent cells. Although the match between experiment and calculation is not perfect, we can see that the main peaks are well reproduced. Similar agreement may be found when comparing with the experimental result of Finney et al. [7] Differences in internal molecular distribution between experimental and theoretical samples in terms of density are commented on in the next section. They may be responsible for the lack of complete agreement. 3.4. Porosity Porosity is one of the most important physical parameters for astrophysical ices [3,31,32]. It is defined by the formula
p=1−
ρaverage ρintrinsic
(1)
where ρaverage is the average density of the ice (including pores and cavities in its structure) and ρintrinsic its intrinsic density (excluding pores). Experimentally, it is difficult to determine porosity from this formula, because direct measurements of both densities require special procedures. Fig. 2. Energy and density variation upon optimization with free cell parameters, for Tamorph = 40 K (above) and Tamorph = 140 K (below). Squares represent (Eaverage, ρ) of the fixed-cell optimized structures and circles represent the corresponding values with released cell parameters. Zero energy for the drawings is set at −5564.4 eV.
3.4.1. Average density Amorphous ices, when formed experimentally, contain voids, pores and other defects. Ices formed by, e.g. vapor-deposition are well known to contain pores. Their measured density is therefore what we call average density. There is a limited number of laboratories with the adequate equipment to measure this quantity, which includes a Quartz Crystal Microbalance (QCMB), to weigh the mass deposited per unit area, and an interferometric technique to measure the ice thickness [33,34], as used e.g. by Satorre et al. [35] Average density can also be determined by means of alternative techniques [30,36,37], and indirectly, from the porosity using the Lorentz-Lorenz relationship [38–40], although this procedure may lead to inaccurate low density values [41]. The actual results depend on several experimental factors, like temperature, pressure, type of vapor deposition [13], but they usually agree within the corresponding error bars.
amorphization and optimization processes at fixed volume, also drops from initial values of tenths of GPa to values lower than 1 kPa. On the other hand, density change is clearly dependent on the initial density value. From ρ = 0.6 g cm−3 to ρ ≤ 1.15 g cm−3, the cell volume diminishes to increase the density; from ρ > 1.20 g cm−3 onwards the volume expands to reduce the density. This variation is practically linear along the series of calculations, with zero variation at ρ ~ 1.17 g cm−3, close to the experimental value for HAD [28,29], or that obtained by Mishima et al. [9] after relaxing a pressurized ice that had reached a density of 1.31 g cm−3. Again our model seems to be in good agreement with experiments.
3.4.2. Intrinsic density The intrinsic density can be deduced from the Pair Correlation Function of the molecules in the solid structure, obtained by diffraction experiments (neutron, X-ray, electrons), easier to apply to thick perfect crystalline ices.41 This procedure limits the studies to certain crystalline structures, for example to hexagonal ice in the case of water. As far as we know, there are only a few other molecules for which these data on crystalline structures are available. However, the amorphous ices represented by our samples are created by optimization of randomly generated structures, which implies the absence of pores in their interior. In this sense, the densities that we use in our approach can be considered intrinsic densities. The possibility of extending the methodology developed in this paper to other
3.3. Pair Correlation Functions (PCF) Information on the relative positions of molecules and atoms is comprised in the theoretical structure through the PCF, frequently denoted as g(r). We focus on OeO pairs to compare with X-ray diffraction results. Fig. 3 shows in panels 1 and 2 the PCF of two of our samples, those of the lowest energy of the sets with ρ = 1.1 g cm−3 and 1.3 g cm−3 respectively, generated at Tamorph = 40 K. In panel 3 we represent the experimental result of Bizid et al. [19], adapted from Jenniskens et al. [3] The HDA trace corresponds to pressurized ice at 77 K, whereas the LDA trace is for the same ice after warming to 140 K and cooling. 3
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
Fig. 3. PCF of OeO separations. Panels 1 and 2: optimized geometry samples generated at 40 K, of density 1.1 g cm−3 and 1.3 g cm−3 respectively; panel 3: experimental graph, adapted from section b of Fig. 4 of Jenniskens et al. [3]; panel 4: composition of panel 2 and the results of Narten et al. [35] (thin grey). Table 1 Approximate OeO distances in Å measured from experimental (Ref. [3]) and theoretical (this work) PCF graphs. Exp., LDA
Calc., 1.1 g cm−3
Exp., HDA
Calc., 1.3 g cm−3
Assignment
2.8
2.8
2.8
2.7
– 4.5 7
– 4.0 ~6.7
3.5 4.6 6.7
3.4 4.2 6.7
OeH⋯O (H bonds) Interstitial O Non-bond Adjacent cells
Table 2 Calculated porosity for the ices generated in this work using average densities from the literature.
astrophysical relevant ices confers an added interest to our results. In summary, in most previous works, the porosity of amorphous ices is calculated from their corresponding average density, when available, and the intrinsic density of the crystalline structure as a substitute of actual values for amorphous samples. This implies large uncertainties in the results. In our case, we can substitute in Eq. (1) experimental average densities from the literature, and our own values for the intrinsic densities for the most favorable structures taken from Fig. 1. The results are collected in Table 2. Cazaux et al. [22] quote a porosity estimation for HDA at low temperatures in the 0.30–0.40 range, which agrees well with our value using an intrinsic density of 1.1 g cm−3. Our value of 1.0 g cm−3 density at 140 K is also close to previous results. On the other hand, our intrinsic density of 1.3 g cm−3, leads to 0.32 porosity, which is large compared to other authors’ estimations.
ρintrinsic/g cm−3 (this work)
ρaverage/g cm−3 (various authors)
Porosity
T = 40 K
1.1 1.3 1.5
0.67a, 0.73b 0.67, 0.73 0.67, 0.73
0.40, 0.34 0.48, 0.44 0.55, 0.51 0.30–0.40c
T = 140 K
1.0 1.3
0.88a, 0.91d 0.88, 0.92
0.12, 0.09 0.32, 0.29 0.13e
a Satorre et al. [35] experimental value obtained for this work by using the experimental method described in Ref. [26]. b Berland et al. [1]. c Cazaux et al. [22], theoretical calculations. d Brown et al. [42]. e Westley et al. [2].
3.5. Formation of the crystal Given the apparent success of this method to study properties of amorphous ices, we can use the outcome of the method to speculate on other interesting features, like the process of building a crystal from the amorphous structure. We can assume that this process will be most favorable when the internal structure of the amorphous ice is closer to that of the crystalline, in terms of interatomic or intermolecular distances, which is tantamount to density. When ices are formed, molecules impinge on a cold surface, and, 4
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
ρ = 1.15 g cm−3 and at ρ = 0.9 g cm−3, correspond to fairly stable amorphous phases, near the experimental values of HDA and LDA. We may mention in passing that the CASTEP program used for all of these calculations is able to detect any increase in symmetry induced during an optimization process, which may take place if the totally amorphous initial structure reaches at least some degree of symmetry. One such case was found at the density of 1.4 g cm−3 for 40 K. This case, and the minimum in Fig. 4 at ~1.45 g cm−3 do not seem to bear any similarity to any experimental observation. However, based on the apparent success of our calculations for other structures, we dare raising a question regarding the existence of some low-temperature, very high density (~1.4 or 1.45 g cm−3) metastable form, which has yet to be detected. 4. Conclusions We have derived a conceptually simple method to study physical properties of amorphous ices. We make quantum calculations at DFT level on samples contained in a box treated as unit cell for solid state methods. The calculations include the generation of amorphous samples, the optimization of their structure, the study of their energetic properties, and the variation of their density upon releasing the size of the cell. With these tools we can harvest information on the stability of the samples vs their density, their porosity and their PCF characteristics. We have applied this study to the water molecule at two temperatures: 40 K, the maximum temperature for the formation of HDA water ice, and 140 K, close to the formation of crystalline cubic water ice Ic. We have followed the evolution of the energetic properties of amorphous ices with density, over a wide range of density values, from 0.6 g cm−3 to 1.6 g cm−3. Our results are in good agreement with experimental observations in general terms. We detect the most favorable amorphous structures at densities ~0.9 g cm−3, 1.15 g cm−3 and ~1.3 g cm−3. Besides these, our calculations predict the possibility of another form of amorphous ice of density ~1.45 g cm−3. The densities generated in this study are intrinsic densities, and thus provide useful data for the calculation of the porosity. Again, the porosities estimated with our results are in good agreement with previous evaluations for water. The fact that such data do not exist for other molecules gives a special interest to the present method, which could be easily applied to species like CH4, NH3, CH3OH or even mixtures of molecules. Furthermore, the optimized structures achieved by this method can be used to predict the IR spectrum of the corresponding sample. Composite or averaged spectra obtained from the samples calculated for each density could be directly compared to the corresponding experimental spectra, where the presence of broad absorption bands is expected for amorphous species.
Fig. 4. Gibbs free energy for the lowest energy structures in each set of (Tamorph, ρ). Squares, Tamorph = 40 K; circles, Tamorph = 140 K. Zero energy for the drawings is set at −5556.55 eV. The dotted lines are drawn only to guide the eye.
depending on their kinetic energy, they may be immediately trapped where they land, or use some excess energy to find some other more favorable environment. Thus, depending on the experimental conditions, ices will be formed of different density, or indeed they will form isles or clusters of a range of densities. If the molecules have a certain amount of mobility they will naturally tend to segregate in the most stable domains or regions, i.e. those designated by minima in the energy vs density representations. Even more, with large enough mobility for a given temperature, they would tend to reorganize in more stable, crystalline structures, through an exothermic process. We can search within our samples for some special cases which may hint to an especially favorable crystallization process. Thus, we have looked for the most stable structure within each set of ρ and Tamorph, and, by making use of the calculated vibrational frequencies of those structures, we estimate their Gibbs free energies. They are collected in Fig. 4, labelled as Glow, and tabulated in Table S2 of the Supplementary Material. We have found that, in a number of cases, some structures have an energy clearly below that of the average value of the set of ten. The pattern of values for 140 K is close to that of the Eaverage values in Fig. 1, i.e. with two minima at 1.0 and 1.3 g cm−3. On the other hand, the minima for Tamorph = 40 K present several interesting features. The deepest minimum appears at ρ = 1.15 g cm−3, again almost exactly the experimental value for HDA, and there are also two shallow minima at 0.9 g cm−3, not far from the experimental result for LDA, and at ~1.45 g cm−3 higher than any experimental measurement that we know of. We can thus infer from Fig. 4 and the above paragraphs that the most stable amorphous structures at 140 K, of density ~1.0 g cm−3 and ~1.3 g cm−3, possess the most favorable molecular distribution. At a low density of 1.0 g cm−3, the molecules are sufficiently apart from each other to rearrange for the crystallization process with the mobility afforded at 140 K. Since the density of crystalline cubic ice at 140 K has been measured as 0.94 g cm−3, we can expect that water molecules in a domain of 0.95–1.0 g cm−3 density would naturally tend to arrange into the crystalline structure Ic. This implies a drop in the Gibbs free energy of ~0.55 eV (the Gibbs free energy calculated for Ic is −5558.6677 eV). On the other hand, at 1.3 g cm−3, the molecules may be already too close to permit mobility, and a crystalline state may not be reached. This structure may be related to that of the “restrained ice”. [43] The mobility in the solid at 40 K must be very limited, and thus any possibility for the molecules to rearrange into different structures is scant. Therefore, the main minima for that temperature in Fig. 4, at
Credit author statement All three authors have contributed equally in all the relevant aspects of the project. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The calculations have been carried out at Trueno Computing facilities (CSIC). Funds from the Spanish MINECO/FEDER FIS2016-77726C3-1-P and C3-3-P, and MINECO/FEDER CTQ2015-65033-P projects 5
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
are acknowledged. We are grateful to B. Maté and C. Santonja for helpful discussions. [22]
Appendix A. Supplementary material [23]
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.cplett.2020.137222.
[24]
References [25]
[1] B.S. Berland, D.E. Brown, M.A. Tolbert, S.M. George, Refractive index and density of vapor-deposited ice, Geophys. Res. Lett. 22 (24) (1995) 3493–3496, https://doi. org/10.1038/nn.4285. [2] M.S. Westley, G.A. Baratta, R.A. Baragiola, Density and index of refraction of water ice films vapor deposited at low temperatures, J. Chem. Phys. 108 (8) (1998) 3321–3326. [3] P. Jenniskens, D.F. Blake, M.A. Wilson, A. Pohorille, High-density amorphous ice, the frost on interstellar grains, ApJ 455 (1995) 389–401. [4] C.A. Tulk, C.J. Benmore, J. Urquidi, D.D. Klug, J. Neuefeind, B. Tomberli, P.A. Egelstaff, Structural studies of several distinct metastable forms of amorphous ice, Science (80-.) 297 (August) (2002) 1320–1324. [5] T. Loerting, K. Winkel, M. Seidl, M. Bauer, C. Mitterdorfer, P.H. Handle, C.G. Salzmann, E. Mayer, J.L. Finney, D.T. Bowron, How many amorphous ices are there? Phys. Chem. Chem. Phys. 13 (19) (2011) 8783–8794, https://doi.org/10. 1039/c0cp02600j. [6] K. Winkel, D.T. Bowron, T. Loerting, E. Mayer, J.L. Finney, Relaxation effects in low density amorphous ice: two distinct structural states observed by neutron diffraction, J. Chem. Phys. 130 (20) (2009) 204502, https://doi.org/10.1063/1.3139007. [7] J.L. Finney, D.T. Bowron, A.K. Soper, T. Loerting, E. Mayer, A. Hallbrucker, Structure of a new dense amorphous ice, Phys. Rev. Lett. 89 (20) (2002) 11–14, https://doi.org/10.1103/PhysRevLett. 89.205503. [8] O. Mishima, Relationship between melting and amorphization of ice, Nature 384 (1996) 546–549. [9] O. Mishima, L.D. Calvert, E. Whalley, “Melting Ice” I at 77 K and 10 Kbar: a new method of making amorphous solids, Lett. Nat. 310 (1984) 393–395. [10] R.J. Nelmes, J.S. Loveday, T. Strässle, C.L. Bull, M. Guthrie, G. Hamel, S. Klotz, Annealed high-density amorphous ice under pressure, Nat. Phys. 2 (6) (2006) 414–418, https://doi.org/10.1038/nphys313. [11] D.T. Limmer, D. Chandler, Theory of amorphous ices, Proc. Natl. Acad. Sci. U. S. A. 111 (26) (2014) 9413–9418, https://doi.org/10.1073/pnas.1407277111. [12] A. Kouchi, T. Hama, Y. Kimura, H. Hidaka, R. Escribano, N. Watanabe, Matrix sublimation method for the formation of high-density amorphous ice, Chem. Phys. Lett. 658 (2016) 287–292, https://doi.org/10.1016/j.cplett.2016.06.066. [13] Z. Dohnálek, G.A. Kimmel, P. Ayotte, R.S. Smith, B.D. Kay, The deposition angledependent density of amorphous solid water films, J. Chem. Phys. 118 (1) (2003) 364–372, https://doi.org/10.1063/1.1525805. [14] J.L. Finney, A. Hallbrucker, I. Kohl, A.K. Soper, D.T. Bowron, Structures of high and low density amorphous ice by neutron diffraction, Phys. Rev. Lett. 88 (22) (2002) 225503/1–225503/4, https://doi.org/10.1103/PhysRevLett.88.225503. [15] B. Guillot, Y. Guissani, Investigation of vapor-deposited amorphous ice and irradiated ice by molecular dynamics simulation, J. Chem. Phys. 120 (9) (2004) 4366–4382, https://doi.org/10.1063/1.1644095. [16] P.C. Gómez, R. Escribano, Vibrational spectra and physico-chemical properties of astrophysical analogs, Phys. Chem. Chem. Phys. 19 (39) (2017) 26582–26588, https://doi.org/10.1039/c7cp04695b. [17] R. Escribano, P.C. Gómez, B. Maté, G. Molpeceres, E. Artacho, Prediction of the near-IR spectra of ices by: ab initio molecular dynamics, Phys. Chem. Chem. Phys. 21 (18) (2019) 9433–9440, https://doi.org/10.1039/c9cp00857h. [18] P. Miró, C.J. Cramer, Water clusters to nanodrops: a tight-binding density functional study, Phys. Chem. Chem. Phys. 15 (6) (2013) 1837–1843, https://doi.org/ 10.1039/c2cp43305b. [19] A. Bizid, L. Bosio, A. Defrain, M. Oumezzine, Structure of high-density amorphous water. I. X-ray diffraction study, J. Chem. Phys. 87 (4) (1987) 2225–2230, https:// doi.org/10.1063/1.453149. [20] A.H. Narten, C.G. Venkatesh, S.A. Rice, Diffraction pattern and structure of amorphous solid water at 10 and 77 °K, J. Chem. Phys. 64 (3) (1976) 1106–1121, https://doi.org/10.1063/1.432298. [21] B.R. Brooks, C. Brooks, A.D. Mackerell, L. Nilsson, R.J. Petrella, B. Roux, Y. Won,
[26]
[27]
[28]
[29]
[30] [31]
[32]
[33]
[34]
[35]
[36] [37]
[38]
[39]
[40] [41]
[42]
[43]
6
G. Archontis, C. Bartels, S. Boresch, et al., CHARMM: molecular dynamics simulation package, J. Comput. Chem. 30 (10) (2009) 1545–1614, https://doi.org/10. 1002/jcc.21287.CHARMM. S. Cazaux, J.-B. Bossa, H. Linnartz, A.G.G.M. Tielens, Pore evolution in interstellar ice analogues, Astron. Astrophys. 573 (2015) A16, https://doi.org/10.1051/00046361/201424466. M. Praprotnik, D. Janežič, Molecular dynamics integration and molecular vibrational theory. III. The infrared spectrum of water, J. Chem. Phys. 122 (17) (2005), https://doi.org/10.1063/1.1884609. S. Yoo, E. Aprà, X.C. Zeng, S.S. Xantheas, High-level ab initio electronic structure calculations of water clusters (H2O)16 and (H2O)17: a new global minimum for (H2O)16, J. Phys. Chem. Lett. 1 (20) (2010) 3122–3127, https://doi.org/10.1021/ jz101245s. A. Lenz, L. Ojamäe, Computational studies of the stability of the (H2O)100 nanodrop, J. Mol. Struct. THEOCHEM 944 (1–3) (2010) 163–167, https://doi.org/10. 1016/j.theochem.2009.12.033. S.J. Clark, M.D. Segall, C.J. Pickard, P.J. Hasnip, M.I.J. Probert, K. Refson, M.C. Payne, First principles methods using CASTEP, Zeitschrift für Krist. 220 (2005) 567–570, https://doi.org/10.1524/zkri.220.5.567.65075. P. Jenniskens, D. Blake, Structural transitions in amorphous water ice and astrophysical implications, Science (80-.) 265 (5173) (1994) 753–756, https://doi.org/ 10.1126/science.11539186. H. Yurimoto, L. Piani, S. Tachibana, K. Murata, Y. Oba, Y. Endo, N. Watanabe, I. Sugawara, Y. Kimura, T. Hama, et al., Liquid-like behavior of UV-irradiated interstellar ice analog at low temperatures, Sci. Adv. 3 (9) (2017), https://doi.org/10. 1126/sciadv.aao2538 eaao2538. C.G. Venkatesh, S.A. Rice, J.B. Bates, A Raman spectral study of amorphous solid water, J. Chem. Phys. 63 (3) (1975) 1065–1071, https://doi.org/10.1063/1. 431447. J.A. Ghormley, C.J. Hochanadel, Amorphous ice: density and reflectivity, Science 171 (3966) (1971) 62–64, https://doi.org/10.1126/science.171.3966.62. J. Hessinger, B.E. White, R.O. Pohl, Elastic properties of amorphous and crystalline ice films, Planet. Space Sci. 44 (9) (1996) 937–944, https://doi.org/10.1016/00320633(95)00112-3. E.H. Mitchell, U. Raut, B.D. Teolis, R.A. Baragiola, Porosity effects on crystallization kinetics of amorphous solid water: implications for cold icy objects in the outer solar system, Icarus 285 (2017) 291–299, https://doi.org/10.1016/j.icarus.2016. 11.004. M.J. Loeffler, M.H. Moore, P.A. Gerakines, The effects of experimental conditions on the refractive index and density of low-temperature ices: solid carbon dioxide, Astrophys. J. 827 (2) (2016) 98, https://doi.org/10.3847/0004-637X/827/2/98. B.A. Seiber, B.E. Wood, A.M. Smith, P.R. Müller, Density of low temperature ice, Science (80-.) 170 (3958) (1970) 652–654, https://doi.org/10.1126/science.170. 3958.652. M.Á. Satorre, R. Luna, C. Millán, M. Domingo, C. Santonja, Density of ices of astrophysical interest, Laboratory Astrophysics, 2018, pp. 51–69, , https://doi.org/ 10.1007/978-3-319-90020-9_4. H.G. Roe, W.M. Grundy, Buoyancy of ice in the CH4–N2 system, Icarus 219 (2012) 733–736, https://doi.org/10.1016/j.icarus.2012.04.007. M. Bouilloud, N. Fray, Y. Benilan, H. Cottin, M.-C. Gazeau, A. Jolly, Bibliographic review and new measurements of the infrared band strengths of pure molecules at 25 K: H2O, CO2, CO, CH4, NH3, CH3OH, HCOOH and H2CO, Mon. Not. R. Astron. Soc. 451 (2) (2015) 2145–2160, https://doi.org/10.1093/mnras/stv1021. D. Fulvio, B. Sivaraman, G.A. Baratta, M.E. Palumbo, N.J. Mason, Novel measurements of refractive index, density and mid-infrared integrated band strengths for solid O2, N2O and NO2:N2O4 mixtures, Spectrochim. Acta – Part A Mol. Biomol. Spectrosc. 72 (5) (2009) 1007–1013, https://doi.org/10.1016/j.saa.2008.12.030. R. Brunetto, M.A. Barucci, E. Dotto, G. Strazzulla, Ion irradiation of frozen methanol, methane, and benzene: linking to the colors of centaurs and trans-neptunian objects, Astrophys. J. 644 (2006) 646–650, https://doi.org/10.1086/503359. B.S. Berland, D.R. Haynes, M.A. Tolbert, S.M. George, Refractive indices of amorphous and crystalline 203 (1994) 4358–4364. C. Millán, C. Santonja, M. Domingo, R. Luna, M.Á. Satorre, An experimental test for effective medium approximations (EMAs), Astron. Astrophys. 628 (2019) A63, https://doi.org/10.1051/0004-6361/201935153. D.E. Brown, S.M. George, C. Huang, E.K.L. Wong, K.B. Rider, R.S. Smith, B.D. Kay, H2O condensation coefficient and refractive index for vapor-deposited ice from molecular beam and optical interference measurements, J. Phys. Chem. 100 (12) (1996) 4988–4995, https://doi.org/10.1021/jp952547j. P. Jenniskens, D.F. Blake, Crystallization of amorphous water ice in the solar system, ApJ 473 (1996) 1104–1113.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Research paper
Density and porosity of amorphous water ice by DFT methods a,⁎
b
Pedro C. Gómez , Miguel Á. Satorre Aznar , Rafael Escribano a b c
T
c
Departamento de Química Física, Facultad de CC. Químicas, Universidad Complutense, and Unidad Asociada Physical Chemistry UCM-CSIC, 28040 Madrid, Spain Departamento de Física Aplicada, EPS Alcoi, Centro de Tecnologías Físicas, Universitat Politècnica de València, Placeta Ferrándiz-Carbonell s/n, 03801 Alcoi, Spain Instituto de Estructura de la Materia, IEM-CSIC, and Unidad Asociada Physical Chemistry UCM-CSIC, Serrano 123, 28006 Madrid, Spain
H I GH L IG H T S
water ice studied at Density Functional Theory level. • Amorphous (LDA) and high-density (HDA) water ices studied at 40 K and 140 K. • Low-density with experimental results on density, porosity and internal structure. • Agreement • Possibility to apply this method to study i.e. astrophysical ices of other species.
A R T I C LE I N FO
A B S T R A C T
Keywords: Amorphous ices Density Porosity DFT methods
A theoretical method is developed to study amorphous ices, based on solid state calculations at DFT level. The method is applied to amorphous water ice of a large range of densities, between 0.6 g cm−3 and 1.6 g cm−3 to represent low-density and high-density amorphous ices, at 40 K and 140 K. The most stable forms appear for densities in good agreement with experimental data, ρ ~ 0.9–1.0 g cm−3 and ρ ~ 1.1–1.2 g cm−3. The method allows evaluating the porosity through the intrinsic density, providing a possible way to estimate the porosity of astrophysical ices.
1. Introduction
2. Methodology
There may be already dozens of high value scientific papers dealing with amorphous water ice (usually denominated ASW) and its physical properties. These include, to quote just a few, density [1,2], which comprises low density, high density and very high density amorphous ice (LDA, HDA and VHDA, respectively); internal structure including Hbonding [3–8]; energetics and stability [4,9–11]; experimental production [9,12,13]; or theoretical models [14–18]. Within the latter topic, the only justification for one more paper may arise if that paper conveys some new methodology allowing the estimation of not sufficiently understood properties of the ice. Our aim in this investigation is thus to present a high-level theoretical methodology at moderate computational cost, test its validity on some well-known examples, and apply it to study ice porosity, which is clearly an underdeveloped subject, especially for astrophysical ices. We can also speculate on the application of this method to deal with an extremely interesting issue, that of the formation of a crystal.
The approach used in this investigation is based in the consideration of an amorphous solid as a 3-dimensional repetitive unit cell where a number of molecules are located in an unsystematic way. The disorderly arrangement is achieved by unlocking the symmetry of a crystalline structure through a fast force-field molecular dynamics using the Amorphous Cell modulus of the Materials Studio (MS) package. In our models, we use a cubic unit cell containing 12 water molecules. The density of the sample is thus determined by the length of the cell side a. This amorphization process has been carried out at two temperatures, Tamorph = 40 K and 140 K, and at a broad range of density values, for ρ between 0.6 g cm−3 and 1.6 g cm−3 (cubic cell side 8.426 Å and 6.073 Å, respectively); this range of densities spans a wider interval than normally achieved experimentally. Bearing in mind that the most abundant OeO distances in amorphous samples have been measured at ~2.7 and 4.3 Å [19,20], the size of our samples guarantees in every case that the most important molecular interactions are enclosed. The size of the cell and the number of molecules within, are very important factors. We have chosen 12 molecules inside the cell, which
⁎
Corresponding author. E-mail address: [email protected] (P.C. Gómez).
https://doi.org/10.1016/j.cplett.2020.137222 Received 21 December 2019; Received in revised form 10 February 2020; Accepted 12 February 2020 Available online 15 February 2020 0009-2614/ © 2020 Elsevier B.V. All rights reserved.
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
is not a large number but still allows describing intermolecular interactions, like H-bonding and van der Waals forces, using a fairly high level of quantum chemistry theory. Samples with hundreds of water molecules have been studied by much less demanding methods, like force-field dynamics models [21–23]. Quantum chemistry models of higher theoretical level are sometimes used for larger number of molecules, 32, 64 or even higher [24,25], but they demand the use of large computation facilities. With our model, we detect the formation of strong H-bonds among two or three molecules within a cell, which often continue in a chain-like form with molecules in the adjacent cells for higher density cases. In our samples, we have observed water molecules with usual coordination numbers from two to four, but including coordination number five as well in some cases [7]. These facts support the choice of 12 molecules per cell to be large enough to describe the main interactions within the ice. Moreover, we can apply our method to big amount of samples to cover an ample range of densities. The total number of samples that we have considered amounts to 230 at fixed cell geometry, and the same number for free cell calculations, discussed below. Our models do not aim to provide an accurate replica of the experiments, as indeed no model can be, but allow studying interesting physical properties. The validity of this methodology to predict new information must be tested by comparing to available experimental results, as we are doing in this investigation. Geometrical structures are optimized looking for a minimum in their potential energy surface (PES). This process is done at Density Functional Theory (DFT) level by means of the CASTEP [26] modulus of the MS package. CASTEP uses pseudopotentials and plane-wave basis sets. The k-point is modified along the increasing accuracy of our calculations. It ranges from (1 × 1 × 1), with k-point separation of 0.1 Å−1 at the initial stages to (2 × 2 × 2), with separation of 0.07 Å−1 at the final optimization process. We have chosen the Generalized Gradient Approximation (GGA) functional by Perdew-Burke-Ernzerhof (PBE) including Grimme’s D2 dispersion correction, which we have found to be very important when dealing with water molecules [16], and a plane-wave basis set cutoff of 830 eV.
Fig. 1. Electronic energy vs density for the set of selected densities. Each energy corresponds to the average of 10 optimized structures generated at 40 K (squares) and 140 K (circles). Zero energy for the drawing is set at −5564.4 eV, with 140 K values offset by −1.0 eV. The dotted line is only drawn to guide the eye.
The graph for 140 K has two local low values, at ~1.0 and ~1.3 g cm−3. There is some likelihood between both graphs, but the distribution of energies around the minima is different, and the high density local minimum for T = 40 K at ρ ~ 1.5 cm−3, is missing at T = 140 K. In the latter case, the structures in the 0.9–1.1 g cm−3 range are energetically close, whereas the local low at 1.2–1.3 g cm−3 is more marked. The interest of these results arises from the assumed inverse relationship between electronic energy and stability of the sample. Thus, we can expect the most stable amorphous structures to be formed at densities close to the low values shown in Fig. 1. Amorphous experimental ices are formed by isles or regions of water aggregates of different densities which give an overall total density of the sample. We can assume that the densities that look more favorable in our calculations will be dominant within the island structure of the experimental ice. The experimental measurements give densities of 1.1–1.2 g cm−3 for HDA grown at T < 40 K [3,28,29], which can be as high as 1.31 g cm−3 at ice-pressurized experiments.[9] On the other hand, for higher temperature experiments, between 70 K and up to near 130 K, LDA is generated with density close to that of cubic ice, 0.94 g cm−3 [20,30]. We can see that our calculations predict the most favorable structures to be formed with densities near the experimental observations, within the expected uncertainties.
3. Results 3.1. Average energy vs density For each Tamorph and ρ we generate 10 amorphous samples which are relaxed to optimize their structure. We choose 10 samples to ensure diversity of initial structures, and after the optimization, we average their energy and use this result, Eaverage, as a representation of the Tamorph vs ρ set. Fig. 1 collects the corresponding results. A full table with the energy and density of each relaxed structure is given in Table S1 of the Supplementary Material. The temperature of Tamorph = 40 K was chosen because it is considered to be the maximum temperature for the formation of HDA water ice.[27] In turn, Tamorph = 140 K was selected because it is close to the formation of crystalline cubic water ice Ic. These two temperatures are sufficiently far apart to permit a better appreciation of differences in the results.
3.2. Cell release It is also interesting to check the effects of releasing the size and shape of the cell together with the molecular arrangement. In this process, the cell will shrink or expand, depending on the internal stress of the structure, searching for a lower energy minimum, and yielding a new volume and therefore a new density. Thus, we have obtained another set of data, which now includes free-cell energies, Efree-cell, and released densities, ρfree-cell. The results are shown in Fig. 2 for Tamorph = 40 K (top) and 140 K (bottom). As above, average values over 10 structures are plotted in every case. Arrows join (Eaverage, ρ) of fixed-cell optimized structures, as in the previous section (squares), with (Efree-cell, ρfree-cell) of free-cell reoptimized structures (circles). The energy always decreases, sometimes as much as 1 eV, and the pressure, understood as the internal stress created in the cell during the
Tamorph = 40 K Fig. 1 indicates that there exist some values of the density which are energetically more favorable than their neighbours. At 40 K, they are located at ρ ~ 1.1 g cm−3, ~1.3 g cm−3, and ρ ~ 1.5 cm−3. The last one is clearly less promising than those in the lower density region. The Eaverage values in the 1.2–1.4 g cm−3 range are fairly close (see also Table S1 for the precise values), indicating that structures with these density values are almost isoenergetic, with the one at ρ ~ 1.1 g cm−3 being only 0.06 eV higher (1.4 kcal mol−1).
Tamorph = 140 K 2
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
Finally, in panel 4, we present a composition of our PCF for ρ = 1.3 g cm−3 and the trace obtained from the values tabulated by Narten et al. [20] To facilitate the comparison, Table 1 collects the measured values of OeO distances from each trace, and their assignment in terms of OeO characteristics. The most prominent peak, i.e. the most frequently occurring OeO distance, always corresponds to the separation between two O atoms linked via H-bond, O⋯HeO. The peak at ~4.5 Å, interpreted as OeO distance between non-bound OeO pairs, is stronger in less-crowded LDA samples, and weakens at higher densities, where water molecules are closer to each other, developing a peak at ~3.5 Å, ascribed to an interstitial O atom linked to other molecules by weak dispersion forces, giving a 5-coordination number O atom [14,19,20]. Finally, the last outstanding peak, at ~6.7 Å, corresponds to distant OeO pairs, which in our calculations are found between atoms in adjacent cells. Although the match between experiment and calculation is not perfect, we can see that the main peaks are well reproduced. Similar agreement may be found when comparing with the experimental result of Finney et al. [7] Differences in internal molecular distribution between experimental and theoretical samples in terms of density are commented on in the next section. They may be responsible for the lack of complete agreement. 3.4. Porosity Porosity is one of the most important physical parameters for astrophysical ices [3,31,32]. It is defined by the formula
p=1−
ρaverage ρintrinsic
(1)
where ρaverage is the average density of the ice (including pores and cavities in its structure) and ρintrinsic its intrinsic density (excluding pores). Experimentally, it is difficult to determine porosity from this formula, because direct measurements of both densities require special procedures. Fig. 2. Energy and density variation upon optimization with free cell parameters, for Tamorph = 40 K (above) and Tamorph = 140 K (below). Squares represent (Eaverage, ρ) of the fixed-cell optimized structures and circles represent the corresponding values with released cell parameters. Zero energy for the drawings is set at −5564.4 eV.
3.4.1. Average density Amorphous ices, when formed experimentally, contain voids, pores and other defects. Ices formed by, e.g. vapor-deposition are well known to contain pores. Their measured density is therefore what we call average density. There is a limited number of laboratories with the adequate equipment to measure this quantity, which includes a Quartz Crystal Microbalance (QCMB), to weigh the mass deposited per unit area, and an interferometric technique to measure the ice thickness [33,34], as used e.g. by Satorre et al. [35] Average density can also be determined by means of alternative techniques [30,36,37], and indirectly, from the porosity using the Lorentz-Lorenz relationship [38–40], although this procedure may lead to inaccurate low density values [41]. The actual results depend on several experimental factors, like temperature, pressure, type of vapor deposition [13], but they usually agree within the corresponding error bars.
amorphization and optimization processes at fixed volume, also drops from initial values of tenths of GPa to values lower than 1 kPa. On the other hand, density change is clearly dependent on the initial density value. From ρ = 0.6 g cm−3 to ρ ≤ 1.15 g cm−3, the cell volume diminishes to increase the density; from ρ > 1.20 g cm−3 onwards the volume expands to reduce the density. This variation is practically linear along the series of calculations, with zero variation at ρ ~ 1.17 g cm−3, close to the experimental value for HAD [28,29], or that obtained by Mishima et al. [9] after relaxing a pressurized ice that had reached a density of 1.31 g cm−3. Again our model seems to be in good agreement with experiments.
3.4.2. Intrinsic density The intrinsic density can be deduced from the Pair Correlation Function of the molecules in the solid structure, obtained by diffraction experiments (neutron, X-ray, electrons), easier to apply to thick perfect crystalline ices.41 This procedure limits the studies to certain crystalline structures, for example to hexagonal ice in the case of water. As far as we know, there are only a few other molecules for which these data on crystalline structures are available. However, the amorphous ices represented by our samples are created by optimization of randomly generated structures, which implies the absence of pores in their interior. In this sense, the densities that we use in our approach can be considered intrinsic densities. The possibility of extending the methodology developed in this paper to other
3.3. Pair Correlation Functions (PCF) Information on the relative positions of molecules and atoms is comprised in the theoretical structure through the PCF, frequently denoted as g(r). We focus on OeO pairs to compare with X-ray diffraction results. Fig. 3 shows in panels 1 and 2 the PCF of two of our samples, those of the lowest energy of the sets with ρ = 1.1 g cm−3 and 1.3 g cm−3 respectively, generated at Tamorph = 40 K. In panel 3 we represent the experimental result of Bizid et al. [19], adapted from Jenniskens et al. [3] The HDA trace corresponds to pressurized ice at 77 K, whereas the LDA trace is for the same ice after warming to 140 K and cooling. 3
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
Fig. 3. PCF of OeO separations. Panels 1 and 2: optimized geometry samples generated at 40 K, of density 1.1 g cm−3 and 1.3 g cm−3 respectively; panel 3: experimental graph, adapted from section b of Fig. 4 of Jenniskens et al. [3]; panel 4: composition of panel 2 and the results of Narten et al. [35] (thin grey). Table 1 Approximate OeO distances in Å measured from experimental (Ref. [3]) and theoretical (this work) PCF graphs. Exp., LDA
Calc., 1.1 g cm−3
Exp., HDA
Calc., 1.3 g cm−3
Assignment
2.8
2.8
2.8
2.7
– 4.5 7
– 4.0 ~6.7
3.5 4.6 6.7
3.4 4.2 6.7
OeH⋯O (H bonds) Interstitial O Non-bond Adjacent cells
Table 2 Calculated porosity for the ices generated in this work using average densities from the literature.
astrophysical relevant ices confers an added interest to our results. In summary, in most previous works, the porosity of amorphous ices is calculated from their corresponding average density, when available, and the intrinsic density of the crystalline structure as a substitute of actual values for amorphous samples. This implies large uncertainties in the results. In our case, we can substitute in Eq. (1) experimental average densities from the literature, and our own values for the intrinsic densities for the most favorable structures taken from Fig. 1. The results are collected in Table 2. Cazaux et al. [22] quote a porosity estimation for HDA at low temperatures in the 0.30–0.40 range, which agrees well with our value using an intrinsic density of 1.1 g cm−3. Our value of 1.0 g cm−3 density at 140 K is also close to previous results. On the other hand, our intrinsic density of 1.3 g cm−3, leads to 0.32 porosity, which is large compared to other authors’ estimations.
ρintrinsic/g cm−3 (this work)
ρaverage/g cm−3 (various authors)
Porosity
T = 40 K
1.1 1.3 1.5
0.67a, 0.73b 0.67, 0.73 0.67, 0.73
0.40, 0.34 0.48, 0.44 0.55, 0.51 0.30–0.40c
T = 140 K
1.0 1.3
0.88a, 0.91d 0.88, 0.92
0.12, 0.09 0.32, 0.29 0.13e
a Satorre et al. [35] experimental value obtained for this work by using the experimental method described in Ref. [26]. b Berland et al. [1]. c Cazaux et al. [22], theoretical calculations. d Brown et al. [42]. e Westley et al. [2].
3.5. Formation of the crystal Given the apparent success of this method to study properties of amorphous ices, we can use the outcome of the method to speculate on other interesting features, like the process of building a crystal from the amorphous structure. We can assume that this process will be most favorable when the internal structure of the amorphous ice is closer to that of the crystalline, in terms of interatomic or intermolecular distances, which is tantamount to density. When ices are formed, molecules impinge on a cold surface, and, 4
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
ρ = 1.15 g cm−3 and at ρ = 0.9 g cm−3, correspond to fairly stable amorphous phases, near the experimental values of HDA and LDA. We may mention in passing that the CASTEP program used for all of these calculations is able to detect any increase in symmetry induced during an optimization process, which may take place if the totally amorphous initial structure reaches at least some degree of symmetry. One such case was found at the density of 1.4 g cm−3 for 40 K. This case, and the minimum in Fig. 4 at ~1.45 g cm−3 do not seem to bear any similarity to any experimental observation. However, based on the apparent success of our calculations for other structures, we dare raising a question regarding the existence of some low-temperature, very high density (~1.4 or 1.45 g cm−3) metastable form, which has yet to be detected. 4. Conclusions We have derived a conceptually simple method to study physical properties of amorphous ices. We make quantum calculations at DFT level on samples contained in a box treated as unit cell for solid state methods. The calculations include the generation of amorphous samples, the optimization of their structure, the study of their energetic properties, and the variation of their density upon releasing the size of the cell. With these tools we can harvest information on the stability of the samples vs their density, their porosity and their PCF characteristics. We have applied this study to the water molecule at two temperatures: 40 K, the maximum temperature for the formation of HDA water ice, and 140 K, close to the formation of crystalline cubic water ice Ic. We have followed the evolution of the energetic properties of amorphous ices with density, over a wide range of density values, from 0.6 g cm−3 to 1.6 g cm−3. Our results are in good agreement with experimental observations in general terms. We detect the most favorable amorphous structures at densities ~0.9 g cm−3, 1.15 g cm−3 and ~1.3 g cm−3. Besides these, our calculations predict the possibility of another form of amorphous ice of density ~1.45 g cm−3. The densities generated in this study are intrinsic densities, and thus provide useful data for the calculation of the porosity. Again, the porosities estimated with our results are in good agreement with previous evaluations for water. The fact that such data do not exist for other molecules gives a special interest to the present method, which could be easily applied to species like CH4, NH3, CH3OH or even mixtures of molecules. Furthermore, the optimized structures achieved by this method can be used to predict the IR spectrum of the corresponding sample. Composite or averaged spectra obtained from the samples calculated for each density could be directly compared to the corresponding experimental spectra, where the presence of broad absorption bands is expected for amorphous species.
Fig. 4. Gibbs free energy for the lowest energy structures in each set of (Tamorph, ρ). Squares, Tamorph = 40 K; circles, Tamorph = 140 K. Zero energy for the drawings is set at −5556.55 eV. The dotted lines are drawn only to guide the eye.
depending on their kinetic energy, they may be immediately trapped where they land, or use some excess energy to find some other more favorable environment. Thus, depending on the experimental conditions, ices will be formed of different density, or indeed they will form isles or clusters of a range of densities. If the molecules have a certain amount of mobility they will naturally tend to segregate in the most stable domains or regions, i.e. those designated by minima in the energy vs density representations. Even more, with large enough mobility for a given temperature, they would tend to reorganize in more stable, crystalline structures, through an exothermic process. We can search within our samples for some special cases which may hint to an especially favorable crystallization process. Thus, we have looked for the most stable structure within each set of ρ and Tamorph, and, by making use of the calculated vibrational frequencies of those structures, we estimate their Gibbs free energies. They are collected in Fig. 4, labelled as Glow, and tabulated in Table S2 of the Supplementary Material. We have found that, in a number of cases, some structures have an energy clearly below that of the average value of the set of ten. The pattern of values for 140 K is close to that of the Eaverage values in Fig. 1, i.e. with two minima at 1.0 and 1.3 g cm−3. On the other hand, the minima for Tamorph = 40 K present several interesting features. The deepest minimum appears at ρ = 1.15 g cm−3, again almost exactly the experimental value for HDA, and there are also two shallow minima at 0.9 g cm−3, not far from the experimental result for LDA, and at ~1.45 g cm−3 higher than any experimental measurement that we know of. We can thus infer from Fig. 4 and the above paragraphs that the most stable amorphous structures at 140 K, of density ~1.0 g cm−3 and ~1.3 g cm−3, possess the most favorable molecular distribution. At a low density of 1.0 g cm−3, the molecules are sufficiently apart from each other to rearrange for the crystallization process with the mobility afforded at 140 K. Since the density of crystalline cubic ice at 140 K has been measured as 0.94 g cm−3, we can expect that water molecules in a domain of 0.95–1.0 g cm−3 density would naturally tend to arrange into the crystalline structure Ic. This implies a drop in the Gibbs free energy of ~0.55 eV (the Gibbs free energy calculated for Ic is −5558.6677 eV). On the other hand, at 1.3 g cm−3, the molecules may be already too close to permit mobility, and a crystalline state may not be reached. This structure may be related to that of the “restrained ice”. [43] The mobility in the solid at 40 K must be very limited, and thus any possibility for the molecules to rearrange into different structures is scant. Therefore, the main minima for that temperature in Fig. 4, at
Credit author statement All three authors have contributed equally in all the relevant aspects of the project. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The calculations have been carried out at Trueno Computing facilities (CSIC). Funds from the Spanish MINECO/FEDER FIS2016-77726C3-1-P and C3-3-P, and MINECO/FEDER CTQ2015-65033-P projects 5
Chemical Physics Letters 745 (2020) 137222
P.C. Gómez, et al.
are acknowledged. We are grateful to B. Maté and C. Santonja for helpful discussions. [22]
Appendix A. Supplementary material [23]
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.cplett.2020.137222.
[24]
References [25]
[1] B.S. Berland, D.E. Brown, M.A. Tolbert, S.M. George, Refractive index and density of vapor-deposited ice, Geophys. Res. Lett. 22 (24) (1995) 3493–3496, https://doi. org/10.1038/nn.4285. [2] M.S. Westley, G.A. Baratta, R.A. Baragiola, Density and index of refraction of water ice films vapor deposited at low temperatures, J. Chem. Phys. 108 (8) (1998) 3321–3326. [3] P. Jenniskens, D.F. Blake, M.A. Wilson, A. Pohorille, High-density amorphous ice, the frost on interstellar grains, ApJ 455 (1995) 389–401. [4] C.A. Tulk, C.J. Benmore, J. Urquidi, D.D. Klug, J. Neuefeind, B. Tomberli, P.A. Egelstaff, Structural studies of several distinct metastable forms of amorphous ice, Science (80-.) 297 (August) (2002) 1320–1324. [5] T. Loerting, K. Winkel, M. Seidl, M. Bauer, C. Mitterdorfer, P.H. Handle, C.G. Salzmann, E. Mayer, J.L. Finney, D.T. Bowron, How many amorphous ices are there? Phys. Chem. Chem. Phys. 13 (19) (2011) 8783–8794, https://doi.org/10. 1039/c0cp02600j. [6] K. Winkel, D.T. Bowron, T. Loerting, E. Mayer, J.L. Finney, Relaxation effects in low density amorphous ice: two distinct structural states observed by neutron diffraction, J. Chem. Phys. 130 (20) (2009) 204502, https://doi.org/10.1063/1.3139007. [7] J.L. Finney, D.T. Bowron, A.K. Soper, T. Loerting, E. Mayer, A. Hallbrucker, Structure of a new dense amorphous ice, Phys. Rev. Lett. 89 (20) (2002) 11–14, https://doi.org/10.1103/PhysRevLett. 89.205503. [8] O. Mishima, Relationship between melting and amorphization of ice, Nature 384 (1996) 546–549. [9] O. Mishima, L.D. Calvert, E. Whalley, “Melting Ice” I at 77 K and 10 Kbar: a new method of making amorphous solids, Lett. Nat. 310 (1984) 393–395. [10] R.J. Nelmes, J.S. Loveday, T. Strässle, C.L. Bull, M. Guthrie, G. Hamel, S. Klotz, Annealed high-density amorphous ice under pressure, Nat. Phys. 2 (6) (2006) 414–418, https://doi.org/10.1038/nphys313. [11] D.T. Limmer, D. Chandler, Theory of amorphous ices, Proc. Natl. Acad. Sci. U. S. A. 111 (26) (2014) 9413–9418, https://doi.org/10.1073/pnas.1407277111. [12] A. Kouchi, T. Hama, Y. Kimura, H. Hidaka, R. Escribano, N. Watanabe, Matrix sublimation method for the formation of high-density amorphous ice, Chem. Phys. Lett. 658 (2016) 287–292, https://doi.org/10.1016/j.cplett.2016.06.066. [13] Z. Dohnálek, G.A. Kimmel, P. Ayotte, R.S. Smith, B.D. Kay, The deposition angledependent density of amorphous solid water films, J. Chem. Phys. 118 (1) (2003) 364–372, https://doi.org/10.1063/1.1525805. [14] J.L. Finney, A. Hallbrucker, I. Kohl, A.K. Soper, D.T. Bowron, Structures of high and low density amorphous ice by neutron diffraction, Phys. Rev. Lett. 88 (22) (2002) 225503/1–225503/4, https://doi.org/10.1103/PhysRevLett.88.225503. [15] B. Guillot, Y. Guissani, Investigation of vapor-deposited amorphous ice and irradiated ice by molecular dynamics simulation, J. Chem. Phys. 120 (9) (2004) 4366–4382, https://doi.org/10.1063/1.1644095. [16] P.C. Gómez, R. Escribano, Vibrational spectra and physico-chemical properties of astrophysical analogs, Phys. Chem. Chem. Phys. 19 (39) (2017) 26582–26588, https://doi.org/10.1039/c7cp04695b. [17] R. Escribano, P.C. Gómez, B. Maté, G. Molpeceres, E. Artacho, Prediction of the near-IR spectra of ices by: ab initio molecular dynamics, Phys. Chem. Chem. Phys. 21 (18) (2019) 9433–9440, https://doi.org/10.1039/c9cp00857h. [18] P. Miró, C.J. Cramer, Water clusters to nanodrops: a tight-binding density functional study, Phys. Chem. Chem. Phys. 15 (6) (2013) 1837–1843, https://doi.org/ 10.1039/c2cp43305b. [19] A. Bizid, L. Bosio, A. Defrain, M. Oumezzine, Structure of high-density amorphous water. I. X-ray diffraction study, J. Chem. Phys. 87 (4) (1987) 2225–2230, https:// doi.org/10.1063/1.453149. [20] A.H. Narten, C.G. Venkatesh, S.A. Rice, Diffraction pattern and structure of amorphous solid water at 10 and 77 °K, J. Chem. Phys. 64 (3) (1976) 1106–1121, https://doi.org/10.1063/1.432298. [21] B.R. Brooks, C. Brooks, A.D. Mackerell, L. Nilsson, R.J. Petrella, B. Roux, Y. Won,
[26]
[27]
[28]
[29]
[30] [31]
[32]
[33]
[34]
[35]
[36] [37]
[38]
[39]
[40] [41]
[42]
[43]
6
G. Archontis, C. Bartels, S. Boresch, et al., CHARMM: molecular dynamics simulation package, J. Comput. Chem. 30 (10) (2009) 1545–1614, https://doi.org/10. 1002/jcc.21287.CHARMM. S. Cazaux, J.-B. Bossa, H. Linnartz, A.G.G.M. Tielens, Pore evolution in interstellar ice analogues, Astron. Astrophys. 573 (2015) A16, https://doi.org/10.1051/00046361/201424466. M. Praprotnik, D. Janežič, Molecular dynamics integration and molecular vibrational theory. III. The infrared spectrum of water, J. Chem. Phys. 122 (17) (2005), https://doi.org/10.1063/1.1884609. S. Yoo, E. Aprà, X.C. Zeng, S.S. Xantheas, High-level ab initio electronic structure calculations of water clusters (H2O)16 and (H2O)17: a new global minimum for (H2O)16, J. Phys. Chem. Lett. 1 (20) (2010) 3122–3127, https://doi.org/10.1021/ jz101245s. A. Lenz, L. Ojamäe, Computational studies of the stability of the (H2O)100 nanodrop, J. Mol. Struct. THEOCHEM 944 (1–3) (2010) 163–167, https://doi.org/10. 1016/j.theochem.2009.12.033. S.J. Clark, M.D. Segall, C.J. Pickard, P.J. Hasnip, M.I.J. Probert, K. Refson, M.C. Payne, First principles methods using CASTEP, Zeitschrift für Krist. 220 (2005) 567–570, https://doi.org/10.1524/zkri.220.5.567.65075. P. Jenniskens, D. Blake, Structural transitions in amorphous water ice and astrophysical implications, Science (80-.) 265 (5173) (1994) 753–756, https://doi.org/ 10.1126/science.11539186. H. Yurimoto, L. Piani, S. Tachibana, K. Murata, Y. Oba, Y. Endo, N. Watanabe, I. Sugawara, Y. Kimura, T. Hama, et al., Liquid-like behavior of UV-irradiated interstellar ice analog at low temperatures, Sci. Adv. 3 (9) (2017), https://doi.org/10. 1126/sciadv.aao2538 eaao2538. C.G. Venkatesh, S.A. Rice, J.B. Bates, A Raman spectral study of amorphous solid water, J. Chem. Phys. 63 (3) (1975) 1065–1071, https://doi.org/10.1063/1. 431447. J.A. Ghormley, C.J. Hochanadel, Amorphous ice: density and reflectivity, Science 171 (3966) (1971) 62–64, https://doi.org/10.1126/science.171.3966.62. J. Hessinger, B.E. White, R.O. Pohl, Elastic properties of amorphous and crystalline ice films, Planet. Space Sci. 44 (9) (1996) 937–944, https://doi.org/10.1016/00320633(95)00112-3. E.H. Mitchell, U. Raut, B.D. Teolis, R.A. Baragiola, Porosity effects on crystallization kinetics of amorphous solid water: implications for cold icy objects in the outer solar system, Icarus 285 (2017) 291–299, https://doi.org/10.1016/j.icarus.2016. 11.004. M.J. Loeffler, M.H. Moore, P.A. Gerakines, The effects of experimental conditions on the refractive index and density of low-temperature ices: solid carbon dioxide, Astrophys. J. 827 (2) (2016) 98, https://doi.org/10.3847/0004-637X/827/2/98. B.A. Seiber, B.E. Wood, A.M. Smith, P.R. Müller, Density of low temperature ice, Science (80-.) 170 (3958) (1970) 652–654, https://doi.org/10.1126/science.170. 3958.652. M.Á. Satorre, R. Luna, C. Millán, M. Domingo, C. Santonja, Density of ices of astrophysical interest, Laboratory Astrophysics, 2018, pp. 51–69, , https://doi.org/ 10.1007/978-3-319-90020-9_4. H.G. Roe, W.M. Grundy, Buoyancy of ice in the CH4–N2 system, Icarus 219 (2012) 733–736, https://doi.org/10.1016/j.icarus.2012.04.007. M. Bouilloud, N. Fray, Y. Benilan, H. Cottin, M.-C. Gazeau, A. Jolly, Bibliographic review and new measurements of the infrared band strengths of pure molecules at 25 K: H2O, CO2, CO, CH4, NH3, CH3OH, HCOOH and H2CO, Mon. Not. R. Astron. Soc. 451 (2) (2015) 2145–2160, https://doi.org/10.1093/mnras/stv1021. D. Fulvio, B. Sivaraman, G.A. Baratta, M.E. Palumbo, N.J. Mason, Novel measurements of refractive index, density and mid-infrared integrated band strengths for solid O2, N2O and NO2:N2O4 mixtures, Spectrochim. Acta – Part A Mol. Biomol. Spectrosc. 72 (5) (2009) 1007–1013, https://doi.org/10.1016/j.saa.2008.12.030. R. Brunetto, M.A. Barucci, E. Dotto, G. Strazzulla, Ion irradiation of frozen methanol, methane, and benzene: linking to the colors of centaurs and trans-neptunian objects, Astrophys. J. 644 (2006) 646–650, https://doi.org/10.1086/503359. B.S. Berland, D.R. Haynes, M.A. Tolbert, S.M. George, Refractive indices of amorphous and crystalline 203 (1994) 4358–4364. C. Millán, C. Santonja, M. Domingo, R. Luna, M.Á. Satorre, An experimental test for effective medium approximations (EMAs), Astron. Astrophys. 628 (2019) A63, https://doi.org/10.1051/0004-6361/201935153. D.E. Brown, S.M. George, C. Huang, E.K.L. Wong, K.B. Rider, R.S. Smith, B.D. Kay, H2O condensation coefficient and refractive index for vapor-deposited ice from molecular beam and optical interference measurements, J. Phys. Chem. 100 (12) (1996) 4988–4995, https://doi.org/10.1021/jp952547j. P. Jenniskens, D.F. Blake, Crystallization of amorphous water ice in the solar system, ApJ 473 (1996) 1104–1113.