- Email: [email protected]
Disorder in ice polymorphs: A Monte Carlo simulation study
Journal of Non-Crystalline Solids 353 (2007) 2698–2707 www.elsevier.com/locate/jnoncrysol
Disorder in ice polymorphs: A Monte Carlo simulation study Albert Barto´k, Andra´s Baranyai
*
Department of Theoretical Chemistry, Eo¨tvo¨s University, 1518 Budapest 112, P.O. Box 32, Hungary Received 27 October 2006; received in revised form 7 May 2007 Available online 5 July 2007
Abstract Two water molecules connected by hydrogen bond in hexagonal ice can have four possible configurations. These configurations are distinguished by the relative orientation of the two molecules and termed for obvious reasons as c-cis, h-cis, c-trans, and h-trans. The occurrence of symmetry permitted dimer orientations is a characteristic feature of each ice polymorph. In the proton-ordered structures the occurrence of orientations is strictly determined, while in the proton-disordered structures it can vary within certain limits. We performed Monte Carlo simulations using the so-called TIP5P-EW, TIP4P-EW and TIP4P-2005 interaction models to study this isomerism for the polymorphs of ice. We found that the variation of energy with the frequency of different dimer orientations in the proton-disordered phases is large enough to influence the results of phase stability studies. Knowing the distributions of dimer orientations of the ice IX–ice III ordered–disordered polymorph pairs, we could estimate the internal energy of ice IX using dimer energies assigned to certain orientations in the disordered phase of ice III. In agreement with experimental evidences at low temperatures the TIP4PEW and TIP4P-2005 potentials predicted lower energy for ice VIII than for ice VII. Ó 2007 Elsevier B.V. All rights reserved. PACS: 05.20.y; 61.66.f; 64.60.Cn Keywords: Crystals; Neutron diffraction/scattering; X-ray diffraction; Microstructure; Modeling and simulation; Monte Carlo simulations; Phases and equilibria; Pressure effects; Structure; Thermal properties; Thermodynamics; Water
1. Introduction Modeling water in its frozen form provides an additional opportunity to understand the behavior of this vital substance. This is indicated by the increasing number of studies in the area [1–8]. In most cases, the variety of structures is accompanied with a relatively small difference in configuration energy. In this respect, it is mandatory to start the calculations from arrangements representing the real phase as closely as possible. The character of the water molecule requires the fulfillment of the so-called ice rule which is a geometrical prescription for the formation of the perfect hydrogen bond network [9]. Nine of the known crystalline phases have disorder in their hydrogen bond structures [9]. In the case of the perfectly tetrahedral *
Corresponding author. Tel.: +36 1 372 2933; fax: +36 1 372 2592. E-mail address: [email protected] (A. Baranyai).
0022-3093/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2007.05.016
phases, a simple geometrical observation shows that the hydrogen bonds can be classified based on the relative orientation of the two molecules in the bond. In Fig. 1, we present the scheme of the four possible arrangements termed as h-cis, c-cis, h-trans, and c-trans. The molecules of the overlapping oxygens are numbered as 1, in the front (donor), and 2, in the background (acceptor). The nonbonding hydrogen of molecule 1, H(1)2 is shown in the front, a thick vertical line connects it to its oxygen atom. The hydrogen atoms of the acceptor molecule are connected by dashed lines. In the scheme of the h-trans arrangement we show the bisector of the H–O–H angle of molecule 2. The angle between the H(1)2–O(1) line and the bisector will be used later on as an unambiguous measure of mutual orientation. The above described isomerism of hydrogen bonds in ice polymorphs has been known for a long time [10]. The proton-ordered and disordered pairs of polymorph structures
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
H(2)2
H(2)1
H(2)2
2699
H(1)2
H(1)2
H(2)1
h-cis H(2)1
H(2)1
H(1)2
H(2)2
H(1)2
H(2)2
c-cis H(2)1 H(2)1
H(2)2 H(1)2
H(2)1
h-trans
H(1)2
H(2)2
H(2)1
H(2)2
H(2)2
c-trans H(1)2
H(1)2
Fig. 1. A schematic diagram of the four possible arrangements of hydrogen-bonded water dimers in the hexagonal phase of ice. In a row we show the front view (left) and the side view (right) of the same dimer. In the front view the oxygen atoms (gray circle) perfectly overlap. The numbers 1 and 2 in parentheses mark atoms of the donor and the acceptor molecules, respectively. On the diagram of the h-cis (side view) and the c-cis (front view) configurations overlapping hydrogen atoms are shifted deliberately for the sake of clarity. On the diagram of the h-trans configuration, we show the angle used later on for defining the relative orientation of the two molecules unambiguously.
such as ice XI and ice Ih, ice IX and ice III, and ice VIII and ice VII indicate that small energy differences of dimer orientations play an important role in phase stability. A detailed study was performed recently by Buch et al. [11] and Rick [12] focusing primarily on the pairs of ice Ih and ice XI. (We adopted the nomenclature of Buch et al. [11] Note that Rick [12] uses different names for the four possible orientations: c-cis – oblique mirror; c-trans – inverse mirror; h-cis – inverse center; h-trans – oblique center.) For these authors, the puzzling question to answer is the experimental evidence that the ferroelectric Cmc21 structure of ice XI contains only 25% of trans bonds in contrast to the 60% trans configurations in the corre-
sponding disordered phase, ice Ih. Still, below 72 K after doping with alkali hydroxides the hexagonal ice undergoes a transition to form ice XI despite the plausible assumption that the trans conformations are energetically more favorable because the interatomic H–H distances are larger [11]. As an orientation for the reader, using the SPC/E model [13], we calculated the relative energies of the four conformations in the hexagonal ice. Using the assignation of Fig. 1, the following distances are identical for all dimer configurations: O(1)–O(2), O(1)–H(2)1, O(1)–H(2)2, O(2)– H(1)1, O(2)–H(1)2, H(1)1–H(2)1, H(1)1–H(2)2. The only distances to be calculated are the H(1)2–H(2)1 and H(1)2– H(2)2 which determine the energy differences of the four
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
arrangements. The results relative to the c-trans pair using ˚ as O–O distance are as follows: h-cis 8.95 kJ/mol; 2.76 A c-cis 7.08 kJ/mol; h-trans 2.54 kJ/mol. The lower configuration energy of the ferroelectric phase is a special effect caused by the long-ranged forces, modeled in simulations by Coulomb interactions. In a detailed study, Rick showed that the NvdE [11] potential predicts lower interaction energy for the ferroelectric structure than for the disordered hexagonal phase [12]. This was not the case for the TIP5P-E [14] (in our notation TIP5P-EW) potential which provided higher energy for the ferroelectric structure [12]. Obviously, the properties of the water model influence the relative energy of the dimers. The position of the partial charges and the H–O–H angle determine the energy differences. The simulation results of Rick and references therein showed that the results depend very heavily on model properties and details of the actual calculations [12]. The problem of relative stability for ice Ih vs. ice XI was investigated recently by Singer et al. [15] using electronic density functional calculations. Their results were, at least qualitatively, in agreement with experimental evidences [9]. These authors tested their method on the transition of the proton-disordered ice VII to the antiferroelectric ice VIII [15]. In Section 2, we present the results of our calculations using the rigid, tetrahedral TIP5P-EW model of Rick [14], (in his notation TIP5P-E) which is a re-parameterization of the TIP5P [16] model to provide a better fit when long-ranged forces are calculated by Ewald summation. We also apply two recently re-parametrized versions of the rigid, planar TIP4P potential [17]: the TIP4P-EW version of Horn et al. [18], and the TIP4P-2005 of Abascal and Vega [19]. The TIP4P-EW provided slightly better results in terms of density and internal energy for both the liquid [18] and the crystalline phases [2] than the original TIP4P. The TIP4P-2005 was developed primarily for the ice phases [19]. The calculations unless stated otherwise report results obtained by the TIP5P-EW potential. This potential produced larger phase stability than the planar potentials. (For instance, ice VII at 2.4 GPa and at 300 K where its structure was determined experimentally [9], melts for all the planar potentials we studied [1,2].) Since this can be the result of the somewhat larger density produced by the TIP5P-EW model, we crosschecked our calculations by using the TIP4P-EW and TIP4P-2005 potentials. In our study, we focus on other forms of ice with emphasis on the ice III vs. ice IX pair. Since the problem of ice Ih vs. ice XI was studied exhaustively [11,12,15], we involve the problem of proton-order in hexagonal ice only for completeness: to introduce the ideas and to check the accuracy of our calculations. We would like to point out at the start that it is difficult for classical rigid interaction models to capture all the subtleties of real ice polymorphs differing only in their proton order. Nevertheless, there are interesting features of these problems understandable by simulation as well.
2. Details and results of calculation 2.1. Hexagonal ice We performed isotherm–isobar Monte Carlo calculations of Rahman-Parrinello type as described in our previous paper [1,2]. To obtain accurate and reproducible results, we considered the simulation cell of the model as the unit cell of the structure. The positions of oxygens were taken from measurement data [9]. Following the method of Buch et al. [11], we created a starting arrangement of the unit cell. Then using the method of Rick and Haymet [20], we generated 105 starting configurations. Configurations with zero net dipole moment were identified and their cis/trans ratio was determined. The configurations with the highest and smallest cis/trans ratio were selected and used as starting cell configurations for standard Monte Carlo runs. We also chose several samples between these extremes. (Restriction of the net dipole to zero makes the interaction energy with the surrounding dielectric medium, in fact, with the rest of the ice phase zero. Simulating molecular reorientations in ice Ih, Rick [12] allowed nonzero net dipole moments of the system with the aim to calculate the dielectric constant.) We performed 105 trial moves per particle in the Monte Carlo process ensuring that half of the attempts were accepted. The method described above was followed for all crystals we studied. Phases of Ih and VII have tetrahedral O–O bonds, departing from the exact value only by thermal motions. The lattice symmetry of hexagonal ice constraints the possible occurrence of the four configurations. The sum of the h-cis and h-trans arrangements must be 75% for geometrical reasons. In the case of ice Ih, using the unit cell of 432 molecules the smallest amount of trans configurations found was 56.3%, while the largest was 64.7%. If the distribution of protons is random substantially larger differences in cis/trans ratios cannot be generated even for larger unit cells. (We checked this by some pilot calculations.) In Fig. 2, we show one particular distribution of halving 1.8 1.6 1.4
frequency of angles (in arbit rary units)
2700
1.2 1.0 0.8 0.6 0.4 0.2 0.0 0
20
40
60
80
100
120
140
160
180
angle (in degrees) Fig. 2. Frequency of dimer angle distribution for the hexagonal ice at 250 K. Number of particles is 432.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
angles as defined by Fig. 1 for hexagonal ice at 250 K. There are four distinct peaks of the distribution reflecting the h-cis (28.4%), c-cis (16.4%), h-trans (46.6%) and c-trans (8.6%) conformations in this order with increasing angles. (For comparison these values for the fully random lattice are 25%, 16.67%, 50%, and 8.33%, respectively [12].) Clearly, without constraining the net dipole moment to zero very different distributions can be generated. Such is the Cmc2l structure of ice IX with only 25% trans, or the possible antiferroelectric Pna21 structure with 100% trans, or the structure found to produce the lowest energy by Buch et al. [11] with 75% trans configurations. In Fig. 3, we present the energies of the hexagonal phase at 250 K in terms of the trans configurations. (Qualitatively identical behavior was found for these systems at 5 K with much smaller error bars.) The larger the ratio of trans configurations the lower the energy of the hexagonal phase. Knowing the results of Buch et al. [11] and Rick [12] this is not a surprise. However, we do not know what the actual
59.32
Energy(KJ/mol)
59.28 59.24
2701
cis/trans ratio is in a real hexagonal crystal (assuming that this ratio cannot vary a lot in nature). The presently available interaction models are not sophisticated enough for selecting a correct set of possible dipole orientation arrangements on thermodynamic grounds. Although, if one studies problems related to phase stability it is mandatory to use a set of large unit cells representing the average structure of the disordered phase correctly. (A possible compromise could be to choose a configuration identical with or very close to the fully random h-bond arrangement.) This conclusion might seem trivial, still calculations of phase diagrams or phase stability of ice phases use a single, ad hoc realization of proton disorder. In Table 1, we present a summary of results for the disordered ice phases calculated by the TIP5P-EW model [14]. Results obtained by the TIP4P-EW and TIP4P-2005 models show qualitatively identical behavior. We present the highest and the lowest energies and densities found for each ice phase. The range of energy for every case is larger than the error of the calculation and comparable to the energy difference between proton-ordered and disordered or stable and metastable pairs of polymorphs. (Here, we assume that the inaccuracy of the interaction potential has much less impact on the energy differences.) In addition to the measurement conditions, we always show close to 0 K calculations because only these energies can be compared to the internal energy calculations. See Refs. [1,2] for details.
59.20
2.2. Ice III and ice VI 59.16 59.12 0.54
0.56
0.58
0.6 0.62 percentege of trans
0.64
0.66
0.68
Fig. 3. Calculated absolute values of energy for the hexagonal phase at different cis/trans ratios at 250 K.
The classification of Fig. 1, however, cannot be done if the arrangement of the molecules is not perfectly tetrahedral. A substantially different geometry of coordination makes this classification ambiguous. However, this is true only for the cis/trans identification. The actual dimer orientation angles are correct if the hydrogen bond linear
Table 1 Collection of characteristic data of simulation Type
N
Ih
432
III
324
IV
1024
V
756
VI
640
VII
432
XII
648
5K
Measurement conditions 3
Energy range (kJ/mol)
Density range (g/cm )
T (K)
Energy range (kJ/mol)
Density range (g/cm3)
P (GPa)
59.30 59.41 55.68 56.00 55.49 55.81 54.96 55.09 54.13 54.33 50.05 50.19 54.17 54.34
1.042 1.042 1.258 1.270 1.422 1.429 1.410 1.412 1.523 1.527 1.753 1.753 1.449 1.452
250
51.83 51.98 48.30 48.70 52.65 52.96 48.54 48.73 48.52 48.71 43.13 43.27 45.45 47.02
0.974 0.974 1.195 1.202 1.377 1.385 1.322 1.325 1.445 1.446 1.627 1.630 1.334 1.346
0.00
250 110 234 225 300 260
0.28 0.00 0.50 1.10 2.40 0.50
The calculations have been carried out at 5 K and at the measurement temperatures [9]. The measurement pressure was applied at both temperatures (shown in the last column). N is the number of particles. The error of the energy calculation at 5 K is 0.01 kJ/mol, while at the highest used temperature (250–300 K) it is 0.1–0.2 kJ/mol depending on the particular phase. The error in the density is less than 0.002 g/cm3.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
2702
between near neighbors. The geometrical properties of the water model have no impact on the starting distributions because in the beginning the donor protons are placed on the O–O lines. By fixing the position of the bisector, the H–O–H angle of the rigid molecule is increased or decreased till it reaches the prescribed angle of the water model. For phases without perfect tetrahedral structure, the unequivocal grouping of bonds should be replaced by a looser argumentation based on the distribution of angles. This is the case for ice III, ice IV, ice V, ice VI, and ice XII. In these phases the structure of the oxygen lattice is the dominant measurement information. If partial proton ordering is suspected (ice V [9]) or there are O–O–O bond angles substantially different from the tetrahedral one (ice IV, ice V), the resulting distribution is a characteristic of the oxygen sublattice and no quantitative conclusion can be made from model calculations. Partial ordering has been observed in the case of ice III, too [9,21–23]. However, in this paper, we omitted to study this problem and considered the phase of ice III as a fully proton disordered phase. In Fig. 4, we show two distributions analogous to that presented in Fig. 2. Fig. 4(a) represents the dimer orientations of ice III cooled down gradually to 5 K, while
a
1.6
frequency of angles (in arbitrary units)
1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0
20
40
60
80
100 120
140
160
180
140
160
180
angle (in degrees)
b 0.8
frequency of angles (in arbitrary units)
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0
20
40
60
80
100
120
angle (in degrees) Fig. 4. Frequency of dimer angle distribution for ice III. (a) Vertical lines border the regions separated by us (5 K) and (b) the same as Fig. 4(a) at 250 K.
Fig. 4(b) shows the same system at 250 K simulated for a cell of 324 molecules. It is true for all the low-temperature distributions shown in this paper that they are practically identical with the dimer angle distribution of the input. The vertical lines introduced by us in the low-temperature distribution in Fig. 4(a) define the borders for the four different regions of orientations (30°, 95°, 140°). The situation in the case of ice III is a bit less elegant then for ice Ih because there are more peaks on the distribution. However, deducing some information from the equivalence of site positions in the crystal and observing the low temperature distribution, we could carry out the division. The positions of boundaries are not identical with the boundaries of Ih, but the extent of the four regions is similar. (In the following, for all the proton-disordered ice phases discussed the shown dimer angle distributions will represent one particular configuration.) The relative height of the peaks varies among the numerous input arrangements generated. These differences will be exploited in the following. (The energy and density ranges we could identify for different polymorphs are shown in Table 1.) It is remarkable that while the low temperature structure of ice III contains several different dimer orientations, regions of the high temperature distribution are placed as in the hexagonal phase. At this stage the question is whether this remarkable similarity is an inherent property of the ice III structure or a mere artifact caused by the rudimentary model. In the case of the hexagonal crystal, the angle distribution of hydrogen-bonded dimers shows four well-defined curves with gradually increasing width as the temperature increases. However, in the case of ice III, there is a large number of different dimer angles at low temperatures as we can see in Fig. 4(a). The TIP5P-EW model overestimates the densities. For Ih it is 0.975 g/cm3, while for ice III it is 1.199 g/cm3, instead of the measured 0.920 g/cm3 and 1.165 g/cm3, respectively. If we compare the pair-correlation functions of the simulated ice Ih and ice III in Fig. 5(a) and (b), we can see that the simulated samples are clearly distinct at 250 K. Similar curves have been published by Vega et al. [6] It is interesting to see Fig. 1 of Abascal et al. [7] which presents the phase diagram of two models compared to the experimental curves. The simulated phase boundaries are similar in character to the experimental ones but at different values of state variables. The calculated diagram is shifted substantially towards low temperatures. It might be that the simulated ice III at 250 K is in fact a compressed hexagonal ice, close to the boundary of the two phases. We carried out simulations at 200 K to see is there any systematic change at lower temperatures in the paircorrelation functions and the dimer angle distribution of ice III. We noticed, however, no qualitative change in behavior. The peaks became a bit sharper due to the lower temperature but the features remained identical. This means that the hydrogen bond isomerism of ice III is less pronounced but very similar in character to that of the hexagonal ice at higher temperatures.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707 7
6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
6
frequency of angles (in arbitrary units)
g (r)
a
5 4 3 2 1 0
0
1
2
3
4
5
6
7
8
r (Å)
0
20
40
60
80
100
120
140
160
180
angle (in degrees)
b 4.0
Fig. 6. Frequency of dimer angle distribution for ice IX at 5 K.
3.5 3.0 2.5
g (r)
2703
2.0 1.5 1.0 0.5 0.0 0
1
2
3
4
5
6
7
8
r (Å) Fig. 5. Partial pair-correlation functions of ices Ih (a) and ice III (b) at 250 K. Full curve: H–H; dashed curve: O–H; dotted curve: O–O.
Analogously, to ice Ih we varied the relative number of dimer types in the defined ranges of ice III as shown in Fig. 4(a) using the method of Rick and Haymet [20], i.e., we changed the relative heights of the peaks in Fig. 4(a) and (b). (Maintaining the ice rule and the zero net dipole of the system we rearranged the hydrogen bond network to give slightly different distributions.) We assigned energies as unknown variables to the four ranges. Using several possible distributions with different relative peak heights within the four ranges we could construct an overdetermined system of linear equations for the energies. By solving these equations we could relate the resulting energy values to the four ranges of the distribution shown in Fig. 4(a). The energy values at 5 K following the ranges of angles from 0° to 180° were 56.0, 57.3, 53.7, and 58.5 kJ/mol, respectively. The error of the solution was 1.1 kJ/mol. The density of ice III varies a bit as we vary the relative weights of different dimer angles (Table 1). This is the cause of the large error in the solution of the linear equations. Unlike to the hexagonal phase, frequency of dimer angles in the four defined ranges has no sum rules. A further difference between these phases is the hierarchy of assigned energy values. Loosely speaking the third angle range corresponds to the h-trans orientation, still its energy is the smallest.
In Fig. 6, we show the distribution of angles for ice IX at 5 K. The three peaks in this figure have their corresponding peaks in Fig. 4(a). The densities of the two phases at 5 K are not very different (ice III: on the average 1.27 g/cm3; ice IX: 1.275 g/cm3) Therefore, we can estimate the energy of ice IX if we combine the corresponding energy values of ice III. Using the weights of the corresponding ice IX peaks we obtain 56:0 13 58:5 23 ffi 57:7 kJ/mol. The value obtained directly from the ice IX simulation is 57.95 kJ/ mol. Since the estimated density of the TIP5P-EW model is too high, we carried out simulations for the TIP4P-EW [18] and the TIP4P-2005 [19] models to check the validity of the calculations above. These models predict the density of ice III at 250 K and at 0.28 GPa 1.172 ± 0.003 g/cm3 and 1.161 ± 0.003 g/cm3, respectively. These values are good approximations of the experimental value of 1.165 g/cm3. The energy values for TIP4P-EW at 10 K following the same ranges of angles (0°–30°, 30°–95°, 95°–140°) were 68.5, 38.1, 62.1, and 57.5 kJ/mol, respectively. From these data the calculated energy for ice IX is 61.2 kJ/mol, while the simulated energy is 61.2 kJ/mol. The energy values for TIP4P-2005 at 10 K following again the same ranges of angles from 0° to 180° were 70.4, 38.6, 63.6, and 58.6 kJ/mol, respectively. From these data the calculated energy for ice IX is 62.5 kJ/mol, while the simulated energy is 62.3 kJ/mol. These calculations revealed a remarkable agreement. While energy values are functions of several factors (system size, calculation of long range interactions, etc.), the leading term of the energy difference for a proton ordered–disordered phase pair is the occurrence of certain dimer orientation angles. It is important to mention by passing that the imperfect tetrahedral order of ice III and ice IX has no impact on our calculations because their imperfections are identical in character and cancel out any distortions influencing the results. As for the partial ordering of ice III, the results above would hold, if the increased order appeared on the dimer angle distribution as a linear combination of the
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
2704
dimer angle distribution of ice III and ice IX, i.e., the new figure of the distribution was a linear combination of Figs. 4(a) and 6. In Fig. 7, we show the low temperature distribution of dimer angles for ice VI. This sample contained 270 water particles. (This figure was practically identical for a system of 640 particles.) There are 12 well-defined peaks representing the very special structure of this phase [24]. One can see these dimer angle distributions as some kind of ‘fingerprint’ of the structure. At 225 K the distribution in terms of the position of the peaks becomes similar to that of the hexagonal phase as shown in Fig. 2. We carried out a similar procedure for this phase what we did for the ice III phase. The resulting energies by the TIP5P-EW potential at 5 K in the order of increasing angles are 51.9, 54.4, 51.9, and 57.4 kJ/mol. The same values at 225 K are 45.8, 49.1, 45.8, and 52.2 kJ/mol, respectively. We do not believe that dimer energies deduced from overdetermined linear equations represent the actual situation in a real crystal. However, the calculations presented above demonstrate that the occurrence of certain orientations is the dominant factor to determine the internal energy and the relative stability of the phases. In the proton disordered phases the actual distributions are restricted by the ice rule and the net zero dipole moment, therefore, there exist only a limited range of typical orientations for the crystal structure. The high temperature distributions contain much less peaks. The curve of the high temperature distribution for ice VI (not shown here, for space saving reasons) is very similar to that of ice III. 2.3. Ice VII and ice VIII The small difference in energies between these phases did not support the type of calculations we carried out for the ice III–IX polymorph pairs. We think that the presently available simple interaction models do not deserve the computational effort to reach the necessary accuracy to do this. While presenting the figures of the characteristic
distributions of orientations, we focused on a slightly different but closely related problem whether simple rigid potential models are able to predict the sign of the energy difference between these phases correctly. To ensure equilibrium conditions we created the input of the ice VIII phase by introducing its proton order into the starting configuration of ice VII. The cubic system transformed into a tetragonal one, manifesting the properties of the equilibrated ice VIII model phase. Such a way we had identical number of particles. (Using the potential of TIP5P-EW both phases were somewhat denser than the experimental values.) In Fig. 8, we show the dimer angle distribution for ice VII at 300 K. (0°, 120°) Following the nomenclature of the hexagonal phase, in this structure there are only two orientations, h-cis and h-trans. We could generate random orientations in the range of 63.8–74.2% for the trans conformations. In the ice VIII phase all the orientations belong to the h-cis category, i.e., only the first peak at 0° of Fig. 8 is present. The transition of the proton-disordered ice VII to the proton-ordered ice VIII is considered well characterized by Singer et al. [15]. According to the available experimental evidence [25,26] the transition takes place in the pressure range of 2–12 GPa around 263–273 K. The TIP5P-EW model used by us could not reproduce these results. We obtained 48.7 kJ/mol for the proton-ordered phase at 5 K which is higher than the energy of the ice VII phase at this temperature (see Table 1). We repeated the calculations for the ice VII and ice VIII polymorph pairs using the TIP4P-EW and TIP4P-2005 potentials. The results are given in Table 2. The energy of ice VIII was lower for both potentials, at each state point studied. This means that even these simple planar, nonpolarizable models are able to capture the ice VII to ice VIII phase transition correctly. This could be an additional test of performance of empirical potentials. (It should be noted, however, that these potentials cannot stabilize these two phases at temperatures >250 K at the lower part of the pressure range.)
2.0 1.2
1.8 1.6
frequency of angles (in arbitrary units)
frequency of angles (in arbitrary units)
1.0 0.8 0.6 0.4
1.4 1.2 1.0 0.8 0.6 0.4
0.2
0.2 0.0
0.0 0
20
40
60
80
100 120
140 160
180
angle (in degrees) Fig. 7. Frequency of dimer angle distributions for ice VI at 5 K.
0
20
40
60
80
100 120
140
160
180
angle (in degrees) Fig. 8. Frequency of dimer angle distribution for ice VII at 300 K.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707 Table 2 Results of simulations for the ice VII and ice VIII phases using the TIP4PEW or the TIP4P-2005 potentials
10 GPa 10 K
TIP4P-2005 2.4 GPa 10 GPa
200 K
10 K
10 K
2.4 GPa 200 K
10 K
VII energy 44.97 41.78 51.72 46.53 43.30 53.30 VIII energy 45.11 41.87 51.78 46.69 43.44 53.41 VII density 1.868 1.833 1.659 1.833 1.801 1.632 VIII density 1.869 1.836 1.662 1.833 1.801 1.633
(in arbitrary units)
TIP4P-EW
0.3
frequency of angles
Ice phase
Energies in kJ/mol, densities in g/cm3. Number of particles is 432. Typical error is ±0.04 kJ/mol in energies, and ±0.002 g/cm3 in densities.
0.2
0.1
0.0 0
40
60
80
100
120
140
160
180
Fig. 9. Frequency of dimer angle distribution for ice IV at 110 K. TIP4P-2005
2.4 GPa
10 GPa
2.4 GPa
44.32 45.02 45.11 45.12
51.07 51.70 51.78 51.81
45.82 46.58 46.69 46.69
52.62 53.33 53.41 53.42
Temperature is 10 K for all cases. Error bars are given in Table 2.
We carried out an additional test of number dependence. In Table 3 we show results of calculations for the ice VIII phase. In addition to proving that our calculations did not depend on particle numbers, this test pointed out the importance of the proper system size when small energy differences are studied.
0.4
(in arbitrary units)
10 GPa
frequency of angles
216 384 432 600
TIP4P-EW
20
angle (in degrees)
Table 3 Number dependence of the ice VIII energies N
2705
0.3
0.2
0.1
0.0 0
20
40
60
80
100 120 140 160 180
angle (in degrees) Fig. 10. Frequency of dimer angle distribution for ice V at 234 K.
2.4. Ice IV, ice V, and ice XII These structures are the most different from the regular tetrahedral arrangements. The O–O–O bond angles vary between wide ranges which creates nonlinear h-bonds. Therefore, in this chapter the expression ‘dimer angle distribution’ is used just as an analogy with the previous chapters without accurate relevance to hydrogen positions. We include these fairly complex structures for completeness and for a comparison with a really random arrangement, liquid water. In Fig. 9, we present the angle distribution of the dimers for the ice IV phase at 110 K. (50°, 120°) At low temperature the distribution of orientation angles contains quite a few peaks but they can be grouped into the three peaks shown in the high temperature distribution. Interestingly, the minima of this distribution approximately match the maxima of the hexagonal-type distributions presented in the previous subchapter. As the pressure increases the perfect tetrahedral coordination of ice phases gets more and more distorted. This is manifested in the different character of the distribution of Fig. 9. The dimer angle distribution of ice V is very complex at low temperature (45°, 80°, 130°). At 234 K the shape of the curve in Fig. 10 reminds to that of Fig. 9 but more irregu-
lar. The narrowness of energy ranges shown in Table 1 indicates that one cannot form a small number of welldefined and separable orientation groups of water molecule pairs. This is manifested in the unstructured character of the distribution. The possible doubt of the crystalline nature of ice V can be removed by looking at the dimer angle distribution of liquid water simulated at ambient conditions. Fig. 11 was created by selecting oxygen pairs if they were closer to each other than the first minimum of the O–O partial pair-correlation function. They were hydrogen bonded if any of the four hydrogens were closer to the other oxygen than the first minimum on the O–H partial pair-correlation function. If more than one hydrogen were within this limit then the one with the smallest O–H distance was selected. Fulfilling the requirements above the two water molecules were considered as a hydrogen-bonded dimer and their geometrical features were determined analogously to that of their crystalline counterparts. The presence of hydrogen bonds in liquid water is a strongly established fact. Still, these hydrogen bonds do not require long-living mutual dimer orientations of water molecules. In Fig. 12(a) and (b), the dimer angle distributions of ice XII are shown (40°, 120°). This system, metastable in
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
(in arbitrary units)
frequency of angles
2706 0.6
3. Conclusions
0.5
We carried out a set of Monte Carlo simulations using the TIP5P-EW [14], TIP4P-EW [18], TIP4P-2005 [19] models in order to study the impact of hydrogen bond isomerism on the energies of ice polymorphs. We defined an angle (see Fig. 1) to describe the mutual orientation of two water molecules connected by hydrogen bond. We created distributions of these dimer angles to identify the frequency of orientations characteristic of the structure. In fact, these distributions serve as a kind of ‘fingerprint’ of the phase. This ‘fingerprint’ character is unequivocal in perfectly tetrahedral structures like ice Ih and ice XI, ice VII and ice VIII. For the other structures the cis/trans isomerism cannot be exactly defined because the nonlinear hydrogen bonds. Nevertheless, for more regular structures like ice III and ice VI these dimer orientation angle distributions contain useful information. In particular, for ice IX using the results of ice III simulations we could estimate the internal energy of the ordered phase with excellent accuracy. We divided up the dimer angle distribution of ice III for four ranges between the limits of 0° and 180°. We created several samples (maintaining the ice rule and the net zero dipole momentum) where within the same angle ranges the frequency was different. Then using the calculated energy values of these samples we solved an overdetermined linear equation with the result of assigning an energy value to each angle range. Multiplying the energy value of a certain ranges of dimer angles with their occurence in the ordered phase, we estimated the energy of ice IX. The agreement was very good for all the three potential models studied. The calculations demonstrated that the leading term of the energy differences between crystals of the same polymorph is the number of water dimers with its particular mutual orientation. (Certainly, these energy values of molecule pairs are different from the values obtained from calculations for two gas phase molecules.) There are further effects to influence this picture, but their importance is secondary. It might be interesting to study this phenomenon whether it is a general rule for ordered–disordered polymorph pairs or a mere accident for the ice III–IX system. However, we think that the quality of the present molecular models is not good enough to reach a more convincing conclusion. In the case of the ice VII–ice VIII pair the above study was impossible to carry out with acceptable accuracy because of the small energy difference of these crystals. It was just this small difference of energy which prompted us to test the used potential models. The TIP5P-EW result was in contrast with the accepted evidence that the antiferroelectric phase, ice VIII is the more stable form in the pressure range of 2–12 GPa below the temperature of 263–273 K [15,22,23]. The TIP5P-EW provided higher energies for the ice VIII phase. This can be attributed to the overestimated density values. The two planar models, however, correctly estimated lower interaction energies for the ice VIII phase than for the ice VII phase at low
0.4 0.3 0.2 0.1 0.0 0
20
40
60
80
100 120 140 160 180
angle (in degrees) Fig. 11. Frequency of dimer angle distributions for water at 300 K.
a
0.8
(in arbitrary units)
frequency of angles
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0
20
40
60
80
100 120 140 160 180
angle (in degrees)
b
0.35
(in arbitrary units)
frequency of angles
0.30 0.25 0.20 0.15 0.10 0.05 0.00 0
20
40
60
80
100 120 140 160 180
angle (in degrees) Fig. 12. Frequency of dimer angle distribution for ice XII. (a) At 5 K and (b) at 250 K.
nature, is more ordered then ice IV or ice V. The maxima of Fig. 12(b) occurs at the same angles as happens for the former two phases. The complex structures of the last three phases prohibit any simple model study presented in the previous sections. While they significantly differ from liquid water, their hydrogen bond arrangements are fairly irregular to be studied.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
temperatures. This is in accordance with the reasonably accurate density values predicted by these models. For completeness we created distributions of ‘dimer angles’ for the three less regular ice phases, ice IV, V and XII. Since in these phases a substantial part of the hydrogen bonds is nonlinear, the distributions characterize the O–O sub-lattice indirectly. Then we compared the distributions of these phases to the dimer angle distribution of water. The distribution of water is completely random, so it is different from that of crystalline phases. However, the hydrogen bond structure of ice IV, V, and XII is too complex to reach quantitative conclusion similar to that of ice III–IX or ice VII–VIII. It became clear from this and previous studies [11,12] that the internal energy (and the density) varies in terms of relative orientations of hydrogen-bonded molecule pairs. The size of this variation is not very large (103–10 J/mol), but understandably comparable to energy differences between proton-ordered and proton disordered pairs of polymorphs. This is also the magnitude which determines the energy difference, thus, the relative stability of pairs like ice Ih and Ic. In contrast to present practice of simulation studies, calculating phase diagrams, for instance, cannot be done without using a representative sample of a proton-disordered phase. Acknowledgement A.B. gratefully acknowledges the support of OTKA Grant T043542. References [1] A. Baranyai, A. Barto´k, A.A. Chialvo, J. Chem. Phys. 123 (2005) 54502. [2] A. Baranyai, A. Barto´k, A.A. Chialvo, J. Chem. Phys. 124 (2006) 074507.
2707
[3] E. Sanz, C. Vega, J.L.F. Abascal, L.G. MacDowell, J. Chem. Phys. 121 (2004) 1165. [4] E. Sanz, C. Vega, J.L.F. Abascal, L.G. MacDowell, Phys. Rev. Lett. 92 (2004) 255701. [5] C. McBride, C. Vega, E. Sanz, L.G. MacDowell, J.L.F. Abascal, Mol. Phys. 103 (2005) 1. [6] C. Vega, C. McBride, E. Sanz, J.L.F. Abascal, Phys. Chem. Chem. Phys. 7 (2005) 1450. [7] J.L.F. Abascal, E. Sanz, R. Garcia Fernandez, C. Vega, J. Chem. Phys. 122 (2005) 234511. [8] R.G. Ferna´ndez, J.L.F. Abascal, C. Vega, J. Chem. Phys. 124 (2006) 144506. [9] V.F. Petrenko, R.W. Withworth, Physics of Ice, Oxford University, 1999. [10] N. Bjerrum, Science 115 (1952) 386. [11] V. Buch, P. Sandler, J. Sadlej, J. Phys. Chem. B 102 (1998) 8641. [12] S.W. Rick, J. Chem. Phys. 122 (2005) 94504. [13] H.J. Berendsen, J.R. Grigera, T.P. Straatma, J. Phys. Chem. 91 (1987) 6269. [14] S.W. Rick, J. Chem. Phys. 120 (2004) 6085. [15] S.J. Singer, J.L. Kuo, T.K. Hirsch, C. Knight, L. Ojamae, M.L. Klein, Phys. Rev. Lett. 94 (2005) 135701. [16] M.W. Mahoney, W.L. Jorgensen, J. Chem. Phys. 112 (2000) 8910. [17] W.L. Jorgensen, J. Chandrashekar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [18] H.W. Horn, W.C. Swope, J.W. Pitera, J.D. Madura, T.J. Dick, G.L. Hura, T. Head-Gordon, J. Chem. Phys. 120 (2004) 9665. [19] J.L.F. Abascal, C. Vega, J. Chem. Phys. 123 (2005) 234505. [20] S.W. Rick, A.D.J. Haymet, J. Chem. Phys. 118 (2003) 9291. [21] C. Lobban, J.L. Finney, W.F. Kuhs, J. Chem. Phys. 112 (2000) 7169. [22] J.D. Londono, W.F. Kuhs, J.L. Finney, J. Chem. Phys. 98 (1993) 4878. [23] C. Knight, S.J. Singer, J. Chem. Phys. 125 (2006) 064506. [24] B. Kamb, Science 150 (1965) 205. [25] G.P. Johari, A. Lavergne, E. Whalley, J. Chem. Phys. 61 (1974) 4292; G.P. Johari, A. Lavergne, E. Whalley, J. Chem. Phys. 73 (1980) 4150. [26] W.F. Kuhs, J.L. Finney, C. Vettier, D.V. Bliss, J. Chem. Phys. 81 (1984) 3612.
Disorder in ice polymorphs: A Monte Carlo simulation study Albert Barto´k, Andra´s Baranyai
*
Department of Theoretical Chemistry, Eo¨tvo¨s University, 1518 Budapest 112, P.O. Box 32, Hungary Received 27 October 2006; received in revised form 7 May 2007 Available online 5 July 2007
Abstract Two water molecules connected by hydrogen bond in hexagonal ice can have four possible configurations. These configurations are distinguished by the relative orientation of the two molecules and termed for obvious reasons as c-cis, h-cis, c-trans, and h-trans. The occurrence of symmetry permitted dimer orientations is a characteristic feature of each ice polymorph. In the proton-ordered structures the occurrence of orientations is strictly determined, while in the proton-disordered structures it can vary within certain limits. We performed Monte Carlo simulations using the so-called TIP5P-EW, TIP4P-EW and TIP4P-2005 interaction models to study this isomerism for the polymorphs of ice. We found that the variation of energy with the frequency of different dimer orientations in the proton-disordered phases is large enough to influence the results of phase stability studies. Knowing the distributions of dimer orientations of the ice IX–ice III ordered–disordered polymorph pairs, we could estimate the internal energy of ice IX using dimer energies assigned to certain orientations in the disordered phase of ice III. In agreement with experimental evidences at low temperatures the TIP4PEW and TIP4P-2005 potentials predicted lower energy for ice VIII than for ice VII. Ó 2007 Elsevier B.V. All rights reserved. PACS: 05.20.y; 61.66.f; 64.60.Cn Keywords: Crystals; Neutron diffraction/scattering; X-ray diffraction; Microstructure; Modeling and simulation; Monte Carlo simulations; Phases and equilibria; Pressure effects; Structure; Thermal properties; Thermodynamics; Water
1. Introduction Modeling water in its frozen form provides an additional opportunity to understand the behavior of this vital substance. This is indicated by the increasing number of studies in the area [1–8]. In most cases, the variety of structures is accompanied with a relatively small difference in configuration energy. In this respect, it is mandatory to start the calculations from arrangements representing the real phase as closely as possible. The character of the water molecule requires the fulfillment of the so-called ice rule which is a geometrical prescription for the formation of the perfect hydrogen bond network [9]. Nine of the known crystalline phases have disorder in their hydrogen bond structures [9]. In the case of the perfectly tetrahedral *
Corresponding author. Tel.: +36 1 372 2933; fax: +36 1 372 2592. E-mail address: [email protected] (A. Baranyai).
0022-3093/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2007.05.016
phases, a simple geometrical observation shows that the hydrogen bonds can be classified based on the relative orientation of the two molecules in the bond. In Fig. 1, we present the scheme of the four possible arrangements termed as h-cis, c-cis, h-trans, and c-trans. The molecules of the overlapping oxygens are numbered as 1, in the front (donor), and 2, in the background (acceptor). The nonbonding hydrogen of molecule 1, H(1)2 is shown in the front, a thick vertical line connects it to its oxygen atom. The hydrogen atoms of the acceptor molecule are connected by dashed lines. In the scheme of the h-trans arrangement we show the bisector of the H–O–H angle of molecule 2. The angle between the H(1)2–O(1) line and the bisector will be used later on as an unambiguous measure of mutual orientation. The above described isomerism of hydrogen bonds in ice polymorphs has been known for a long time [10]. The proton-ordered and disordered pairs of polymorph structures
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
H(2)2
H(2)1
H(2)2
2699
H(1)2
H(1)2
H(2)1
h-cis H(2)1
H(2)1
H(1)2
H(2)2
H(1)2
H(2)2
c-cis H(2)1 H(2)1
H(2)2 H(1)2
H(2)1
h-trans
H(1)2
H(2)2
H(2)1
H(2)2
H(2)2
c-trans H(1)2
H(1)2
Fig. 1. A schematic diagram of the four possible arrangements of hydrogen-bonded water dimers in the hexagonal phase of ice. In a row we show the front view (left) and the side view (right) of the same dimer. In the front view the oxygen atoms (gray circle) perfectly overlap. The numbers 1 and 2 in parentheses mark atoms of the donor and the acceptor molecules, respectively. On the diagram of the h-cis (side view) and the c-cis (front view) configurations overlapping hydrogen atoms are shifted deliberately for the sake of clarity. On the diagram of the h-trans configuration, we show the angle used later on for defining the relative orientation of the two molecules unambiguously.
such as ice XI and ice Ih, ice IX and ice III, and ice VIII and ice VII indicate that small energy differences of dimer orientations play an important role in phase stability. A detailed study was performed recently by Buch et al. [11] and Rick [12] focusing primarily on the pairs of ice Ih and ice XI. (We adopted the nomenclature of Buch et al. [11] Note that Rick [12] uses different names for the four possible orientations: c-cis – oblique mirror; c-trans – inverse mirror; h-cis – inverse center; h-trans – oblique center.) For these authors, the puzzling question to answer is the experimental evidence that the ferroelectric Cmc21 structure of ice XI contains only 25% of trans bonds in contrast to the 60% trans configurations in the corre-
sponding disordered phase, ice Ih. Still, below 72 K after doping with alkali hydroxides the hexagonal ice undergoes a transition to form ice XI despite the plausible assumption that the trans conformations are energetically more favorable because the interatomic H–H distances are larger [11]. As an orientation for the reader, using the SPC/E model [13], we calculated the relative energies of the four conformations in the hexagonal ice. Using the assignation of Fig. 1, the following distances are identical for all dimer configurations: O(1)–O(2), O(1)–H(2)1, O(1)–H(2)2, O(2)– H(1)1, O(2)–H(1)2, H(1)1–H(2)1, H(1)1–H(2)2. The only distances to be calculated are the H(1)2–H(2)1 and H(1)2– H(2)2 which determine the energy differences of the four
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
arrangements. The results relative to the c-trans pair using ˚ as O–O distance are as follows: h-cis 8.95 kJ/mol; 2.76 A c-cis 7.08 kJ/mol; h-trans 2.54 kJ/mol. The lower configuration energy of the ferroelectric phase is a special effect caused by the long-ranged forces, modeled in simulations by Coulomb interactions. In a detailed study, Rick showed that the NvdE [11] potential predicts lower interaction energy for the ferroelectric structure than for the disordered hexagonal phase [12]. This was not the case for the TIP5P-E [14] (in our notation TIP5P-EW) potential which provided higher energy for the ferroelectric structure [12]. Obviously, the properties of the water model influence the relative energy of the dimers. The position of the partial charges and the H–O–H angle determine the energy differences. The simulation results of Rick and references therein showed that the results depend very heavily on model properties and details of the actual calculations [12]. The problem of relative stability for ice Ih vs. ice XI was investigated recently by Singer et al. [15] using electronic density functional calculations. Their results were, at least qualitatively, in agreement with experimental evidences [9]. These authors tested their method on the transition of the proton-disordered ice VII to the antiferroelectric ice VIII [15]. In Section 2, we present the results of our calculations using the rigid, tetrahedral TIP5P-EW model of Rick [14], (in his notation TIP5P-E) which is a re-parameterization of the TIP5P [16] model to provide a better fit when long-ranged forces are calculated by Ewald summation. We also apply two recently re-parametrized versions of the rigid, planar TIP4P potential [17]: the TIP4P-EW version of Horn et al. [18], and the TIP4P-2005 of Abascal and Vega [19]. The TIP4P-EW provided slightly better results in terms of density and internal energy for both the liquid [18] and the crystalline phases [2] than the original TIP4P. The TIP4P-2005 was developed primarily for the ice phases [19]. The calculations unless stated otherwise report results obtained by the TIP5P-EW potential. This potential produced larger phase stability than the planar potentials. (For instance, ice VII at 2.4 GPa and at 300 K where its structure was determined experimentally [9], melts for all the planar potentials we studied [1,2].) Since this can be the result of the somewhat larger density produced by the TIP5P-EW model, we crosschecked our calculations by using the TIP4P-EW and TIP4P-2005 potentials. In our study, we focus on other forms of ice with emphasis on the ice III vs. ice IX pair. Since the problem of ice Ih vs. ice XI was studied exhaustively [11,12,15], we involve the problem of proton-order in hexagonal ice only for completeness: to introduce the ideas and to check the accuracy of our calculations. We would like to point out at the start that it is difficult for classical rigid interaction models to capture all the subtleties of real ice polymorphs differing only in their proton order. Nevertheless, there are interesting features of these problems understandable by simulation as well.
2. Details and results of calculation 2.1. Hexagonal ice We performed isotherm–isobar Monte Carlo calculations of Rahman-Parrinello type as described in our previous paper [1,2]. To obtain accurate and reproducible results, we considered the simulation cell of the model as the unit cell of the structure. The positions of oxygens were taken from measurement data [9]. Following the method of Buch et al. [11], we created a starting arrangement of the unit cell. Then using the method of Rick and Haymet [20], we generated 105 starting configurations. Configurations with zero net dipole moment were identified and their cis/trans ratio was determined. The configurations with the highest and smallest cis/trans ratio were selected and used as starting cell configurations for standard Monte Carlo runs. We also chose several samples between these extremes. (Restriction of the net dipole to zero makes the interaction energy with the surrounding dielectric medium, in fact, with the rest of the ice phase zero. Simulating molecular reorientations in ice Ih, Rick [12] allowed nonzero net dipole moments of the system with the aim to calculate the dielectric constant.) We performed 105 trial moves per particle in the Monte Carlo process ensuring that half of the attempts were accepted. The method described above was followed for all crystals we studied. Phases of Ih and VII have tetrahedral O–O bonds, departing from the exact value only by thermal motions. The lattice symmetry of hexagonal ice constraints the possible occurrence of the four configurations. The sum of the h-cis and h-trans arrangements must be 75% for geometrical reasons. In the case of ice Ih, using the unit cell of 432 molecules the smallest amount of trans configurations found was 56.3%, while the largest was 64.7%. If the distribution of protons is random substantially larger differences in cis/trans ratios cannot be generated even for larger unit cells. (We checked this by some pilot calculations.) In Fig. 2, we show one particular distribution of halving 1.8 1.6 1.4
frequency of angles (in arbit rary units)
2700
1.2 1.0 0.8 0.6 0.4 0.2 0.0 0
20
40
60
80
100
120
140
160
180
angle (in degrees) Fig. 2. Frequency of dimer angle distribution for the hexagonal ice at 250 K. Number of particles is 432.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
angles as defined by Fig. 1 for hexagonal ice at 250 K. There are four distinct peaks of the distribution reflecting the h-cis (28.4%), c-cis (16.4%), h-trans (46.6%) and c-trans (8.6%) conformations in this order with increasing angles. (For comparison these values for the fully random lattice are 25%, 16.67%, 50%, and 8.33%, respectively [12].) Clearly, without constraining the net dipole moment to zero very different distributions can be generated. Such is the Cmc2l structure of ice IX with only 25% trans, or the possible antiferroelectric Pna21 structure with 100% trans, or the structure found to produce the lowest energy by Buch et al. [11] with 75% trans configurations. In Fig. 3, we present the energies of the hexagonal phase at 250 K in terms of the trans configurations. (Qualitatively identical behavior was found for these systems at 5 K with much smaller error bars.) The larger the ratio of trans configurations the lower the energy of the hexagonal phase. Knowing the results of Buch et al. [11] and Rick [12] this is not a surprise. However, we do not know what the actual
59.32
Energy(KJ/mol)
59.28 59.24
2701
cis/trans ratio is in a real hexagonal crystal (assuming that this ratio cannot vary a lot in nature). The presently available interaction models are not sophisticated enough for selecting a correct set of possible dipole orientation arrangements on thermodynamic grounds. Although, if one studies problems related to phase stability it is mandatory to use a set of large unit cells representing the average structure of the disordered phase correctly. (A possible compromise could be to choose a configuration identical with or very close to the fully random h-bond arrangement.) This conclusion might seem trivial, still calculations of phase diagrams or phase stability of ice phases use a single, ad hoc realization of proton disorder. In Table 1, we present a summary of results for the disordered ice phases calculated by the TIP5P-EW model [14]. Results obtained by the TIP4P-EW and TIP4P-2005 models show qualitatively identical behavior. We present the highest and the lowest energies and densities found for each ice phase. The range of energy for every case is larger than the error of the calculation and comparable to the energy difference between proton-ordered and disordered or stable and metastable pairs of polymorphs. (Here, we assume that the inaccuracy of the interaction potential has much less impact on the energy differences.) In addition to the measurement conditions, we always show close to 0 K calculations because only these energies can be compared to the internal energy calculations. See Refs. [1,2] for details.
59.20
2.2. Ice III and ice VI 59.16 59.12 0.54
0.56
0.58
0.6 0.62 percentege of trans
0.64
0.66
0.68
Fig. 3. Calculated absolute values of energy for the hexagonal phase at different cis/trans ratios at 250 K.
The classification of Fig. 1, however, cannot be done if the arrangement of the molecules is not perfectly tetrahedral. A substantially different geometry of coordination makes this classification ambiguous. However, this is true only for the cis/trans identification. The actual dimer orientation angles are correct if the hydrogen bond linear
Table 1 Collection of characteristic data of simulation Type
N
Ih
432
III
324
IV
1024
V
756
VI
640
VII
432
XII
648
5K
Measurement conditions 3
Energy range (kJ/mol)
Density range (g/cm )
T (K)
Energy range (kJ/mol)
Density range (g/cm3)
P (GPa)
59.30 59.41 55.68 56.00 55.49 55.81 54.96 55.09 54.13 54.33 50.05 50.19 54.17 54.34
1.042 1.042 1.258 1.270 1.422 1.429 1.410 1.412 1.523 1.527 1.753 1.753 1.449 1.452
250
51.83 51.98 48.30 48.70 52.65 52.96 48.54 48.73 48.52 48.71 43.13 43.27 45.45 47.02
0.974 0.974 1.195 1.202 1.377 1.385 1.322 1.325 1.445 1.446 1.627 1.630 1.334 1.346
0.00
250 110 234 225 300 260
0.28 0.00 0.50 1.10 2.40 0.50
The calculations have been carried out at 5 K and at the measurement temperatures [9]. The measurement pressure was applied at both temperatures (shown in the last column). N is the number of particles. The error of the energy calculation at 5 K is 0.01 kJ/mol, while at the highest used temperature (250–300 K) it is 0.1–0.2 kJ/mol depending on the particular phase. The error in the density is less than 0.002 g/cm3.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
2702
between near neighbors. The geometrical properties of the water model have no impact on the starting distributions because in the beginning the donor protons are placed on the O–O lines. By fixing the position of the bisector, the H–O–H angle of the rigid molecule is increased or decreased till it reaches the prescribed angle of the water model. For phases without perfect tetrahedral structure, the unequivocal grouping of bonds should be replaced by a looser argumentation based on the distribution of angles. This is the case for ice III, ice IV, ice V, ice VI, and ice XII. In these phases the structure of the oxygen lattice is the dominant measurement information. If partial proton ordering is suspected (ice V [9]) or there are O–O–O bond angles substantially different from the tetrahedral one (ice IV, ice V), the resulting distribution is a characteristic of the oxygen sublattice and no quantitative conclusion can be made from model calculations. Partial ordering has been observed in the case of ice III, too [9,21–23]. However, in this paper, we omitted to study this problem and considered the phase of ice III as a fully proton disordered phase. In Fig. 4, we show two distributions analogous to that presented in Fig. 2. Fig. 4(a) represents the dimer orientations of ice III cooled down gradually to 5 K, while
a
1.6
frequency of angles (in arbitrary units)
1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0
20
40
60
80
100 120
140
160
180
140
160
180
angle (in degrees)
b 0.8
frequency of angles (in arbitrary units)
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0
20
40
60
80
100
120
angle (in degrees) Fig. 4. Frequency of dimer angle distribution for ice III. (a) Vertical lines border the regions separated by us (5 K) and (b) the same as Fig. 4(a) at 250 K.
Fig. 4(b) shows the same system at 250 K simulated for a cell of 324 molecules. It is true for all the low-temperature distributions shown in this paper that they are practically identical with the dimer angle distribution of the input. The vertical lines introduced by us in the low-temperature distribution in Fig. 4(a) define the borders for the four different regions of orientations (30°, 95°, 140°). The situation in the case of ice III is a bit less elegant then for ice Ih because there are more peaks on the distribution. However, deducing some information from the equivalence of site positions in the crystal and observing the low temperature distribution, we could carry out the division. The positions of boundaries are not identical with the boundaries of Ih, but the extent of the four regions is similar. (In the following, for all the proton-disordered ice phases discussed the shown dimer angle distributions will represent one particular configuration.) The relative height of the peaks varies among the numerous input arrangements generated. These differences will be exploited in the following. (The energy and density ranges we could identify for different polymorphs are shown in Table 1.) It is remarkable that while the low temperature structure of ice III contains several different dimer orientations, regions of the high temperature distribution are placed as in the hexagonal phase. At this stage the question is whether this remarkable similarity is an inherent property of the ice III structure or a mere artifact caused by the rudimentary model. In the case of the hexagonal crystal, the angle distribution of hydrogen-bonded dimers shows four well-defined curves with gradually increasing width as the temperature increases. However, in the case of ice III, there is a large number of different dimer angles at low temperatures as we can see in Fig. 4(a). The TIP5P-EW model overestimates the densities. For Ih it is 0.975 g/cm3, while for ice III it is 1.199 g/cm3, instead of the measured 0.920 g/cm3 and 1.165 g/cm3, respectively. If we compare the pair-correlation functions of the simulated ice Ih and ice III in Fig. 5(a) and (b), we can see that the simulated samples are clearly distinct at 250 K. Similar curves have been published by Vega et al. [6] It is interesting to see Fig. 1 of Abascal et al. [7] which presents the phase diagram of two models compared to the experimental curves. The simulated phase boundaries are similar in character to the experimental ones but at different values of state variables. The calculated diagram is shifted substantially towards low temperatures. It might be that the simulated ice III at 250 K is in fact a compressed hexagonal ice, close to the boundary of the two phases. We carried out simulations at 200 K to see is there any systematic change at lower temperatures in the paircorrelation functions and the dimer angle distribution of ice III. We noticed, however, no qualitative change in behavior. The peaks became a bit sharper due to the lower temperature but the features remained identical. This means that the hydrogen bond isomerism of ice III is less pronounced but very similar in character to that of the hexagonal ice at higher temperatures.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707 7
6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
6
frequency of angles (in arbitrary units)
g (r)
a
5 4 3 2 1 0
0
1
2
3
4
5
6
7
8
r (Å)
0
20
40
60
80
100
120
140
160
180
angle (in degrees)
b 4.0
Fig. 6. Frequency of dimer angle distribution for ice IX at 5 K.
3.5 3.0 2.5
g (r)
2703
2.0 1.5 1.0 0.5 0.0 0
1
2
3
4
5
6
7
8
r (Å) Fig. 5. Partial pair-correlation functions of ices Ih (a) and ice III (b) at 250 K. Full curve: H–H; dashed curve: O–H; dotted curve: O–O.
Analogously, to ice Ih we varied the relative number of dimer types in the defined ranges of ice III as shown in Fig. 4(a) using the method of Rick and Haymet [20], i.e., we changed the relative heights of the peaks in Fig. 4(a) and (b). (Maintaining the ice rule and the zero net dipole of the system we rearranged the hydrogen bond network to give slightly different distributions.) We assigned energies as unknown variables to the four ranges. Using several possible distributions with different relative peak heights within the four ranges we could construct an overdetermined system of linear equations for the energies. By solving these equations we could relate the resulting energy values to the four ranges of the distribution shown in Fig. 4(a). The energy values at 5 K following the ranges of angles from 0° to 180° were 56.0, 57.3, 53.7, and 58.5 kJ/mol, respectively. The error of the solution was 1.1 kJ/mol. The density of ice III varies a bit as we vary the relative weights of different dimer angles (Table 1). This is the cause of the large error in the solution of the linear equations. Unlike to the hexagonal phase, frequency of dimer angles in the four defined ranges has no sum rules. A further difference between these phases is the hierarchy of assigned energy values. Loosely speaking the third angle range corresponds to the h-trans orientation, still its energy is the smallest.
In Fig. 6, we show the distribution of angles for ice IX at 5 K. The three peaks in this figure have their corresponding peaks in Fig. 4(a). The densities of the two phases at 5 K are not very different (ice III: on the average 1.27 g/cm3; ice IX: 1.275 g/cm3) Therefore, we can estimate the energy of ice IX if we combine the corresponding energy values of ice III. Using the weights of the corresponding ice IX peaks we obtain 56:0 13 58:5 23 ffi 57:7 kJ/mol. The value obtained directly from the ice IX simulation is 57.95 kJ/ mol. Since the estimated density of the TIP5P-EW model is too high, we carried out simulations for the TIP4P-EW [18] and the TIP4P-2005 [19] models to check the validity of the calculations above. These models predict the density of ice III at 250 K and at 0.28 GPa 1.172 ± 0.003 g/cm3 and 1.161 ± 0.003 g/cm3, respectively. These values are good approximations of the experimental value of 1.165 g/cm3. The energy values for TIP4P-EW at 10 K following the same ranges of angles (0°–30°, 30°–95°, 95°–140°) were 68.5, 38.1, 62.1, and 57.5 kJ/mol, respectively. From these data the calculated energy for ice IX is 61.2 kJ/mol, while the simulated energy is 61.2 kJ/mol. The energy values for TIP4P-2005 at 10 K following again the same ranges of angles from 0° to 180° were 70.4, 38.6, 63.6, and 58.6 kJ/mol, respectively. From these data the calculated energy for ice IX is 62.5 kJ/mol, while the simulated energy is 62.3 kJ/mol. These calculations revealed a remarkable agreement. While energy values are functions of several factors (system size, calculation of long range interactions, etc.), the leading term of the energy difference for a proton ordered–disordered phase pair is the occurrence of certain dimer orientation angles. It is important to mention by passing that the imperfect tetrahedral order of ice III and ice IX has no impact on our calculations because their imperfections are identical in character and cancel out any distortions influencing the results. As for the partial ordering of ice III, the results above would hold, if the increased order appeared on the dimer angle distribution as a linear combination of the
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
2704
dimer angle distribution of ice III and ice IX, i.e., the new figure of the distribution was a linear combination of Figs. 4(a) and 6. In Fig. 7, we show the low temperature distribution of dimer angles for ice VI. This sample contained 270 water particles. (This figure was practically identical for a system of 640 particles.) There are 12 well-defined peaks representing the very special structure of this phase [24]. One can see these dimer angle distributions as some kind of ‘fingerprint’ of the structure. At 225 K the distribution in terms of the position of the peaks becomes similar to that of the hexagonal phase as shown in Fig. 2. We carried out a similar procedure for this phase what we did for the ice III phase. The resulting energies by the TIP5P-EW potential at 5 K in the order of increasing angles are 51.9, 54.4, 51.9, and 57.4 kJ/mol. The same values at 225 K are 45.8, 49.1, 45.8, and 52.2 kJ/mol, respectively. We do not believe that dimer energies deduced from overdetermined linear equations represent the actual situation in a real crystal. However, the calculations presented above demonstrate that the occurrence of certain orientations is the dominant factor to determine the internal energy and the relative stability of the phases. In the proton disordered phases the actual distributions are restricted by the ice rule and the net zero dipole moment, therefore, there exist only a limited range of typical orientations for the crystal structure. The high temperature distributions contain much less peaks. The curve of the high temperature distribution for ice VI (not shown here, for space saving reasons) is very similar to that of ice III. 2.3. Ice VII and ice VIII The small difference in energies between these phases did not support the type of calculations we carried out for the ice III–IX polymorph pairs. We think that the presently available simple interaction models do not deserve the computational effort to reach the necessary accuracy to do this. While presenting the figures of the characteristic
distributions of orientations, we focused on a slightly different but closely related problem whether simple rigid potential models are able to predict the sign of the energy difference between these phases correctly. To ensure equilibrium conditions we created the input of the ice VIII phase by introducing its proton order into the starting configuration of ice VII. The cubic system transformed into a tetragonal one, manifesting the properties of the equilibrated ice VIII model phase. Such a way we had identical number of particles. (Using the potential of TIP5P-EW both phases were somewhat denser than the experimental values.) In Fig. 8, we show the dimer angle distribution for ice VII at 300 K. (0°, 120°) Following the nomenclature of the hexagonal phase, in this structure there are only two orientations, h-cis and h-trans. We could generate random orientations in the range of 63.8–74.2% for the trans conformations. In the ice VIII phase all the orientations belong to the h-cis category, i.e., only the first peak at 0° of Fig. 8 is present. The transition of the proton-disordered ice VII to the proton-ordered ice VIII is considered well characterized by Singer et al. [15]. According to the available experimental evidence [25,26] the transition takes place in the pressure range of 2–12 GPa around 263–273 K. The TIP5P-EW model used by us could not reproduce these results. We obtained 48.7 kJ/mol for the proton-ordered phase at 5 K which is higher than the energy of the ice VII phase at this temperature (see Table 1). We repeated the calculations for the ice VII and ice VIII polymorph pairs using the TIP4P-EW and TIP4P-2005 potentials. The results are given in Table 2. The energy of ice VIII was lower for both potentials, at each state point studied. This means that even these simple planar, nonpolarizable models are able to capture the ice VII to ice VIII phase transition correctly. This could be an additional test of performance of empirical potentials. (It should be noted, however, that these potentials cannot stabilize these two phases at temperatures >250 K at the lower part of the pressure range.)
2.0 1.2
1.8 1.6
frequency of angles (in arbitrary units)
frequency of angles (in arbitrary units)
1.0 0.8 0.6 0.4
1.4 1.2 1.0 0.8 0.6 0.4
0.2
0.2 0.0
0.0 0
20
40
60
80
100 120
140 160
180
angle (in degrees) Fig. 7. Frequency of dimer angle distributions for ice VI at 5 K.
0
20
40
60
80
100 120
140
160
180
angle (in degrees) Fig. 8. Frequency of dimer angle distribution for ice VII at 300 K.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707 Table 2 Results of simulations for the ice VII and ice VIII phases using the TIP4PEW or the TIP4P-2005 potentials
10 GPa 10 K
TIP4P-2005 2.4 GPa 10 GPa
200 K
10 K
10 K
2.4 GPa 200 K
10 K
VII energy 44.97 41.78 51.72 46.53 43.30 53.30 VIII energy 45.11 41.87 51.78 46.69 43.44 53.41 VII density 1.868 1.833 1.659 1.833 1.801 1.632 VIII density 1.869 1.836 1.662 1.833 1.801 1.633
(in arbitrary units)
TIP4P-EW
0.3
frequency of angles
Ice phase
Energies in kJ/mol, densities in g/cm3. Number of particles is 432. Typical error is ±0.04 kJ/mol in energies, and ±0.002 g/cm3 in densities.
0.2
0.1
0.0 0
40
60
80
100
120
140
160
180
Fig. 9. Frequency of dimer angle distribution for ice IV at 110 K. TIP4P-2005
2.4 GPa
10 GPa
2.4 GPa
44.32 45.02 45.11 45.12
51.07 51.70 51.78 51.81
45.82 46.58 46.69 46.69
52.62 53.33 53.41 53.42
Temperature is 10 K for all cases. Error bars are given in Table 2.
We carried out an additional test of number dependence. In Table 3 we show results of calculations for the ice VIII phase. In addition to proving that our calculations did not depend on particle numbers, this test pointed out the importance of the proper system size when small energy differences are studied.
0.4
(in arbitrary units)
10 GPa
frequency of angles
216 384 432 600
TIP4P-EW
20
angle (in degrees)
Table 3 Number dependence of the ice VIII energies N
2705
0.3
0.2
0.1
0.0 0
20
40
60
80
100 120 140 160 180
angle (in degrees) Fig. 10. Frequency of dimer angle distribution for ice V at 234 K.
2.4. Ice IV, ice V, and ice XII These structures are the most different from the regular tetrahedral arrangements. The O–O–O bond angles vary between wide ranges which creates nonlinear h-bonds. Therefore, in this chapter the expression ‘dimer angle distribution’ is used just as an analogy with the previous chapters without accurate relevance to hydrogen positions. We include these fairly complex structures for completeness and for a comparison with a really random arrangement, liquid water. In Fig. 9, we present the angle distribution of the dimers for the ice IV phase at 110 K. (50°, 120°) At low temperature the distribution of orientation angles contains quite a few peaks but they can be grouped into the three peaks shown in the high temperature distribution. Interestingly, the minima of this distribution approximately match the maxima of the hexagonal-type distributions presented in the previous subchapter. As the pressure increases the perfect tetrahedral coordination of ice phases gets more and more distorted. This is manifested in the different character of the distribution of Fig. 9. The dimer angle distribution of ice V is very complex at low temperature (45°, 80°, 130°). At 234 K the shape of the curve in Fig. 10 reminds to that of Fig. 9 but more irregu-
lar. The narrowness of energy ranges shown in Table 1 indicates that one cannot form a small number of welldefined and separable orientation groups of water molecule pairs. This is manifested in the unstructured character of the distribution. The possible doubt of the crystalline nature of ice V can be removed by looking at the dimer angle distribution of liquid water simulated at ambient conditions. Fig. 11 was created by selecting oxygen pairs if they were closer to each other than the first minimum of the O–O partial pair-correlation function. They were hydrogen bonded if any of the four hydrogens were closer to the other oxygen than the first minimum on the O–H partial pair-correlation function. If more than one hydrogen were within this limit then the one with the smallest O–H distance was selected. Fulfilling the requirements above the two water molecules were considered as a hydrogen-bonded dimer and their geometrical features were determined analogously to that of their crystalline counterparts. The presence of hydrogen bonds in liquid water is a strongly established fact. Still, these hydrogen bonds do not require long-living mutual dimer orientations of water molecules. In Fig. 12(a) and (b), the dimer angle distributions of ice XII are shown (40°, 120°). This system, metastable in
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
(in arbitrary units)
frequency of angles
2706 0.6
3. Conclusions
0.5
We carried out a set of Monte Carlo simulations using the TIP5P-EW [14], TIP4P-EW [18], TIP4P-2005 [19] models in order to study the impact of hydrogen bond isomerism on the energies of ice polymorphs. We defined an angle (see Fig. 1) to describe the mutual orientation of two water molecules connected by hydrogen bond. We created distributions of these dimer angles to identify the frequency of orientations characteristic of the structure. In fact, these distributions serve as a kind of ‘fingerprint’ of the phase. This ‘fingerprint’ character is unequivocal in perfectly tetrahedral structures like ice Ih and ice XI, ice VII and ice VIII. For the other structures the cis/trans isomerism cannot be exactly defined because the nonlinear hydrogen bonds. Nevertheless, for more regular structures like ice III and ice VI these dimer orientation angle distributions contain useful information. In particular, for ice IX using the results of ice III simulations we could estimate the internal energy of the ordered phase with excellent accuracy. We divided up the dimer angle distribution of ice III for four ranges between the limits of 0° and 180°. We created several samples (maintaining the ice rule and the net zero dipole momentum) where within the same angle ranges the frequency was different. Then using the calculated energy values of these samples we solved an overdetermined linear equation with the result of assigning an energy value to each angle range. Multiplying the energy value of a certain ranges of dimer angles with their occurence in the ordered phase, we estimated the energy of ice IX. The agreement was very good for all the three potential models studied. The calculations demonstrated that the leading term of the energy differences between crystals of the same polymorph is the number of water dimers with its particular mutual orientation. (Certainly, these energy values of molecule pairs are different from the values obtained from calculations for two gas phase molecules.) There are further effects to influence this picture, but their importance is secondary. It might be interesting to study this phenomenon whether it is a general rule for ordered–disordered polymorph pairs or a mere accident for the ice III–IX system. However, we think that the quality of the present molecular models is not good enough to reach a more convincing conclusion. In the case of the ice VII–ice VIII pair the above study was impossible to carry out with acceptable accuracy because of the small energy difference of these crystals. It was just this small difference of energy which prompted us to test the used potential models. The TIP5P-EW result was in contrast with the accepted evidence that the antiferroelectric phase, ice VIII is the more stable form in the pressure range of 2–12 GPa below the temperature of 263–273 K [15,22,23]. The TIP5P-EW provided higher energies for the ice VIII phase. This can be attributed to the overestimated density values. The two planar models, however, correctly estimated lower interaction energies for the ice VIII phase than for the ice VII phase at low
0.4 0.3 0.2 0.1 0.0 0
20
40
60
80
100 120 140 160 180
angle (in degrees) Fig. 11. Frequency of dimer angle distributions for water at 300 K.
a
0.8
(in arbitrary units)
frequency of angles
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0
20
40
60
80
100 120 140 160 180
angle (in degrees)
b
0.35
(in arbitrary units)
frequency of angles
0.30 0.25 0.20 0.15 0.10 0.05 0.00 0
20
40
60
80
100 120 140 160 180
angle (in degrees) Fig. 12. Frequency of dimer angle distribution for ice XII. (a) At 5 K and (b) at 250 K.
nature, is more ordered then ice IV or ice V. The maxima of Fig. 12(b) occurs at the same angles as happens for the former two phases. The complex structures of the last three phases prohibit any simple model study presented in the previous sections. While they significantly differ from liquid water, their hydrogen bond arrangements are fairly irregular to be studied.
A. Barto´k, A. Baranyai / Journal of Non-Crystalline Solids 353 (2007) 2698–2707
temperatures. This is in accordance with the reasonably accurate density values predicted by these models. For completeness we created distributions of ‘dimer angles’ for the three less regular ice phases, ice IV, V and XII. Since in these phases a substantial part of the hydrogen bonds is nonlinear, the distributions characterize the O–O sub-lattice indirectly. Then we compared the distributions of these phases to the dimer angle distribution of water. The distribution of water is completely random, so it is different from that of crystalline phases. However, the hydrogen bond structure of ice IV, V, and XII is too complex to reach quantitative conclusion similar to that of ice III–IX or ice VII–VIII. It became clear from this and previous studies [11,12] that the internal energy (and the density) varies in terms of relative orientations of hydrogen-bonded molecule pairs. The size of this variation is not very large (103–10 J/mol), but understandably comparable to energy differences between proton-ordered and proton disordered pairs of polymorphs. This is also the magnitude which determines the energy difference, thus, the relative stability of pairs like ice Ih and Ic. In contrast to present practice of simulation studies, calculating phase diagrams, for instance, cannot be done without using a representative sample of a proton-disordered phase. Acknowledgement A.B. gratefully acknowledges the support of OTKA Grant T043542. References [1] A. Baranyai, A. Barto´k, A.A. Chialvo, J. Chem. Phys. 123 (2005) 54502. [2] A. Baranyai, A. Barto´k, A.A. Chialvo, J. Chem. Phys. 124 (2006) 074507.
2707
[3] E. Sanz, C. Vega, J.L.F. Abascal, L.G. MacDowell, J. Chem. Phys. 121 (2004) 1165. [4] E. Sanz, C. Vega, J.L.F. Abascal, L.G. MacDowell, Phys. Rev. Lett. 92 (2004) 255701. [5] C. McBride, C. Vega, E. Sanz, L.G. MacDowell, J.L.F. Abascal, Mol. Phys. 103 (2005) 1. [6] C. Vega, C. McBride, E. Sanz, J.L.F. Abascal, Phys. Chem. Chem. Phys. 7 (2005) 1450. [7] J.L.F. Abascal, E. Sanz, R. Garcia Fernandez, C. Vega, J. Chem. Phys. 122 (2005) 234511. [8] R.G. Ferna´ndez, J.L.F. Abascal, C. Vega, J. Chem. Phys. 124 (2006) 144506. [9] V.F. Petrenko, R.W. Withworth, Physics of Ice, Oxford University, 1999. [10] N. Bjerrum, Science 115 (1952) 386. [11] V. Buch, P. Sandler, J. Sadlej, J. Phys. Chem. B 102 (1998) 8641. [12] S.W. Rick, J. Chem. Phys. 122 (2005) 94504. [13] H.J. Berendsen, J.R. Grigera, T.P. Straatma, J. Phys. Chem. 91 (1987) 6269. [14] S.W. Rick, J. Chem. Phys. 120 (2004) 6085. [15] S.J. Singer, J.L. Kuo, T.K. Hirsch, C. Knight, L. Ojamae, M.L. Klein, Phys. Rev. Lett. 94 (2005) 135701. [16] M.W. Mahoney, W.L. Jorgensen, J. Chem. Phys. 112 (2000) 8910. [17] W.L. Jorgensen, J. Chandrashekar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [18] H.W. Horn, W.C. Swope, J.W. Pitera, J.D. Madura, T.J. Dick, G.L. Hura, T. Head-Gordon, J. Chem. Phys. 120 (2004) 9665. [19] J.L.F. Abascal, C. Vega, J. Chem. Phys. 123 (2005) 234505. [20] S.W. Rick, A.D.J. Haymet, J. Chem. Phys. 118 (2003) 9291. [21] C. Lobban, J.L. Finney, W.F. Kuhs, J. Chem. Phys. 112 (2000) 7169. [22] J.D. Londono, W.F. Kuhs, J.L. Finney, J. Chem. Phys. 98 (1993) 4878. [23] C. Knight, S.J. Singer, J. Chem. Phys. 125 (2006) 064506. [24] B. Kamb, Science 150 (1965) 205. [25] G.P. Johari, A. Lavergne, E. Whalley, J. Chem. Phys. 61 (1974) 4292; G.P. Johari, A. Lavergne, E. Whalley, J. Chem. Phys. 73 (1980) 4150. [26] W.F. Kuhs, J.L. Finney, C. Vettier, D.V. Bliss, J. Chem. Phys. 81 (1984) 3612.