Proton ordered cubic and hexagonal periodic models of ordinary ice

Chemical Physics Letters 409 (2005) 110–117 www.elsevier.com/locate/cplett

Proton ordered cubic and hexagonal periodic models of ordinary ice S. Casassa a, M. Calatayud b, K. Doll c, C. Minot

b,*

, C. Pisani

a

a

Dipartimento di Chimica IFM and Centre of Excellence NIS (Nanostructured Interfaces and Surfaces), Universita` di Torino, via Giuria 5, I-10125 Torino, Italy b Laboratoire de Chimie The´orique, Universite´ P. et M. Curie, 4 Place Jussieu case 137, site Ivry, 75252 Paris Cedex 05, France c Institut fu¨r Mathematische Physik, Technische Universita¨t Braunschweig, Mendelssohnstrasse 3, D-38106 Braunschweig, Germany Received 14 February 2005; in final form 19 April 2005 Available online 23 May 2005

Abstract We compare different quantum-mechanical (QM) periodic models of ordinary ice with the objective to provide reference data for use in simulation studies. Four crystalline structures are considered derived from a cubic or hexagonal arrangement of the oxygen sublattice, both polar and apolar, and three Hamiltonians: HF, PW91 and B3LYP. Two periodic programs have been used, CRYSTAL and VASP. For each given QM method, the four structures are almost indistinguishable as far as energy and local bonding structure are concerned. Results are much more dependent on the QM technique concerning relative distances and activation barriers. Ó 2005 Elsevier B.V. All rights reserved.

1. Introduction Water ice is an important object of research in many areas of biological, environmental and atmospheric sciences [1,2], and computer simulations may help us to understand its properties. Therefore, a lot of attention has been recently devoted to the modeling of bulk and surface ice structure and reactivity [3–6]. In many applications, it is mandatory to adopt to this end quantum mechanical (QM) models; this is true in particular when one is interested in investigating the catalytic properties of ice [6–8]. Two main issues determine the adequacy of the model: geometrical structure and level of the QM approximation. Concerning the former issue, cluster models are very frequently adopted, since they permit the use of standard QM molecular packages. There are, however, some inconveniences in this choice. Groups of few water molecules tend to collapse into peculiar structures, very far

*

Corresponding author. E-mail address: [email protected] (C. Minot).

0009-2614/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.04.068

from that of ordinary ice [9]. Furthermore, non-saturated hydrogen bonds at the border may result in strong electrostatic fields inside the cluster, with deep consequences on its chemical properties [10]. As very efficient QM periodic programs are now available, crystalline models of ice can be used instead, which present some attractive features. First of all, a few well defined structures can be characterized once and for all, while an almost infinite variety of cluster models exist. Second, periodic systems analogous to ordinary ice are intrinsically stable; their surfaces, and local defects within them may be modeled and optimized, so reactive properties can be realistically simulated. Nevertheless, ordinary ice simulation by means of periodic models prevents the description of its proton disorder. It is, however, questionable how important disorder may be in determining the chemical properties of ice (local defects are certainly more important). The search for reasonable proton-order is not easily accessible by experiment since spectroscopy mostly characterizes the O positions. More stable structures (proton ordered ice) should appear below a transition

S. Casassa et al. / Chemical Physics Letters 409 (2005) 110–117

temperature. This one has been evaluated to 72 K for KOH-doped ice [11] but 237 K has also been proposed (see [12] and references therein). One should also be careful about the meaning of the apparent disorder. Starting from the local environment of a given O atom, the proton distribution in its surroundings is bound to obey the Bernarl–Fowler rules [13] avoiding any local mismatch. However, there are many ways to distribute hydrogens so as to realize this kind of global ordering, and the different distributions are equally probable. To simulate a fairly weighted set of such possible distributions while preserving translational symmetry would require the use of very large unit cells, and make the model not univocally defined. In any event, only the simplest well-defined proton ordered crystalline structures are considered in the following: the comparison between crystal structures implying unit cells of increasing size (containing 2, 4 and 8 water molecules) can give trends going from the smallest unit cell with the maximal symmetry to larger ones with less symmetry. Proton-order has been previously investigated by theory [14] starting from experimental observations. Ice XI is the only known proton ordered phase of ice Ih [15] system with hexagonal symmetry at zero pressure. The two space groups, CmC21 (mm2) and Pna21 (Pmmm) have been proposed as models for the ice XI. Neutron diffraction measurements [16] and mean field calculations [17] conclude to domains of ordered ferroelectric phase (CmC21). The group space CmC21 has been assigned from accurate experimental measurements at 5 K [18,19]. For the ideal structures, a preference by 0.6 kcal/mol in favor of Pna21 was initially found by ab initio calculations [20]. However, the structures were later on found quasi isoenergetic [14]. From a MonteCarlo simulation, Barkema and de Boer [21] have found a non-ferroelectric ground state (Pna21) and a transition toward a partially ordered phase at T = 36 K, lower than assumed in real ice [16]. Consider now the problem of the QM approach to be adopted. Molecular codes permit the use of a wide variety of Hamiltonians and techniques, ranging from Hartree–Fock (HF) and post-HF, to all kinds of density functional theory (DFT) schemes. Periodic programs are less flexible in this respect, since they allow only one-electron Hamiltonians to be used, namely HF or DFT ones; the so-called ÔhybridÕ approaches are sometimes feasible (see below). As concerns basis sets, molecular programs usually adopt Gaussian type orbitals (GTO) as representative functions; their number and features must be accurately calibrated to obtain the required accuracy. Periodic codes may either use GTOs, or a universal set of plane waves (PW) with the same periodicity of the crystal. In the latter case, the ÔpowerÕ of the basis set is uniquely determined by the PW cutoff that is, their maximum frequency; soft-core pseudopotentials must be used to get rid of the innermost elec-

111

trons, since their description would require an excessive number of PWs. At variance with previous studies, we have concentrated here our attention on two topics: the influence of the crystalline structure, and that of the QM method adopted on the calculated equilibrium properties of the very simplest periodic models of ordinary ice; the objective is to provide reference data which can be useful in modeling.

2. Computational models 2.1. Ice structures All models of ordinary ice must obey the Bernal– Fowler rules [13]: each O atom is at the center of a (distorted) tetrahedron with four other oxygen atoms at its vertices; along all four nearest neighbor O


O direc˚ ) there is an H atom; two of these tions (dOO  2.7 A ˚ ), the two form the OH hydroxyl bond (dO–H  0.95 A ˚ ) with the other O


H hydrogen bond (dO–H  1.75 A central O atom. Once translational symmetry is imposed, equivalent O atoms in different cells have their two OH bonds with the same orientations; since all four tetrahedral directions contain an OH bond, the minimum number of translationally inequivalent water molecules per cell is two. Let us label in this case O 0 , O00 the two inequivalent O species, W 0 , W00 the corresponding water molecules (see Fig. 1). The two space groups where the O 0 and O00 sublattices can interpenetrate in such a way that each O 0 is at the center of a perfect tetrahedron with O00 at its vertices, and vice versa, are

Fig. 1. The minimum number of translationally inequivalent water molecules is two, W 0 (at the center of the figure) and W00 (at the vertices of the tetrahedron shown).

112

S. Casassa et al. / Chemical Physics Letters 409 (2005) 110–117

facepffifficentered cubic (fcc), with lattice parameter ac ¼ ffi ð4= 3Þd OO and hexagonal close packed qffiffi (hcp) with ch = 8dOO/3 and an ideal ratio ah =ch ¼ 38. In the former case, there are two molecules per cell with oxygen atoms in a diamond-like structure; in the latter four, two per orientation. These structures are polar, however since two molecules W 0 and W00 have their electric dipole in the [1 0 0] direction. The same holds true for the hcp case. In order to generate apolar models, the cell must be doubled, at least. Starting from the ideal fcc and hcp structures, in the polar and apolar case, after allowing for tetrahedral distortion and finding the equilibrium configuration, we obtain the four ice structures discussed below: PC, AC, PH, AH (Figs. 2–5). The HOH angle for the ice structure differs from that imposed by the tetrahedral environment so that the optimized structures deviate from the ideal ones. In the following, we will consider ÔidealÕ and ÔoptimizedÕ structures. In the former, the H atoms are forced to be along the axis of a perfect tetrahedron (O 0 O00 directions). Only the OH bond distances (two short and two long) break the tetrahedral symmetry. For all the ÔidealÕ structures, there are only two parameters to optimize: dOO and dd OH . The first one determines the lattice parameters OO and the second one the H positions. The ÔoptimizedÕ structure is fully optimized within the space symmetry group and includes a larger number of parameters. This should allow distinguishing the role of topology or polarity (present in the ideal structure) from that of the degrees of freedom. In the fully ÔoptimizedÕ PC structure, there are four parameters to optimize: at, ct, xH, zH (or at, ct, dd OH and HOH) instead of two. For the other OO

Fig. 3. The apolar cubic lattice.

Fig. 4. The polar hexagonal lattice.

Fig. 2. The polar cubic lattice. The characteristic symmetry element is a screw axis of order 4, parallel to the z axis through the midpoint between nearest neighbor oxygen atoms, and transforming each water molecule into a neighboring one, rotated by 90°.

Fig. 5. The apolar hexagonal lattice.

S. Casassa et al. / Chemical Physics Letters 409 (2005) 110–117

optimized structures, the number of parameters to optimize increases (Table 1) with the loss of symmetry. 2.2. Quantum mechanical techniques Two codes have been used: CRYSTAL [22] and VASP [23]. The CRYSTAL code solves the Schro¨dinger equation for periodic systems both in an HF and in a DFT framework, adopting a set of localized Gaussian-type function. Three Hamiltonians have been tested here: pure HF, pure DFT in the frame of the GGA approximation using the PW91 functional [24], and the ÔhybridÕ B3LYP method, proposed and parameterized by Becke [25], which has met extraordinary success in molecular quantum studies. The basis set adopted is an all-electron 6311G(p,d) one (see Appendix A). With respect to the 6-31G(d) basis set used in previous work on ice structures [10], the present one is appreciably more flexible because of the addition of a polarization function on hydrogens, and of an sp shell (exponent 0.9 a.u.) on oxy-

113

gens. This improvement has been made possible by the new capabilities of the CRYSTAL03 code [22]. Assessment of the quality of the 6311G(p,d) set for the present purposes is provided in Section 2.3, by comparing its performance in describing the water molecule and dimer with that of the simpler 6-31G(d) set, and of a sophisticated 6-311+(2p,2d) set [26] (see Appendix A), which is anyhow too extended for use in periodic calculations. Geometry optimizations have been carried out using the procedure recently implemented by Doll et al. [27]. The binding energy of the four ice models is evaluated with respect to the isolated water molecule; the same basis set of the periodic calculations is used and the geometry optimization is performed for each Hamiltonian (see Table 2). All binding energy data to be reported below, including those for the water dimer, have been corrected for the basis set superposition error (BSSE), by using the a posteriori counterpoise method proposed by Lendvay and Mayer [28]. VASP [23] is a package for performing DFT ab initio quantum-mechanical calculations for periodic systems

Table 1 Space groups symmetry for the ice crystals and relation between the lattice parameters in the ideal structure, referred to ac, the cubic lattice parameter

Hermann-Mauguin symbol Patterson symmetry No. of space group a b c nP

PC

AC

PH

AH

I41md I4/mmm 109 pffiffiffi at ¼ ac =p2ffiffiffi t b ¼ ac = 2 ct = ac 4

P41212 P4/mmm 92 pffiffiffi at ¼ ac =p2ffiffiffi t b ¼ ac = 2 ct = ac 5

Cmc21 (mm2) Cmmm 36 pffiffiffi a ¼ ac p = ffiffiffiffiffiffiffi 2 ffi b ¼ acpffiffiffiffiffiffiffi 3=2ffi c ¼ ac 4=3 12

Pna21 (mm2) Pmmm 33 pffiffiffiffiffiffiffiffi a ¼ ac p4=3 ffiffiffiffiffiffiffiffi b ¼ ac p3=2 ffiffiffi c c¼a = 2 18

nP is the number of independent parameters in the full optimization (for the ideal structure, this number is reduced to 2).

Table 2 Calculated data for the water molecule and the water dimer Hamiltonian

Method

Molecule

Dimer

dOH

ˆH HO

dOH

dOO

dOH/dO–H

Ebind

HF

CRY-a CRY-b CRY-c

0.948 0.941 0.940

105.6 105.4 106.2

0.952 0.945 0.944

2.966 2.983 3.042

0.472 0.464 0.450

4.82 4.36 3.82

B3LYP

CRY-a CRY-b CRY-c

0.969 0.962 0.961

103.7 103.8 105.1

0.977 0.970 0.969

2.872 2.891 2.925

0.513 0.503 0.495

5.29 4.76 4.58

PW91

CRY-a CRY-b CRY-c VASP-a VASP-b

0.976 0.969 0.968 0.973 0.968

103.1 102.9 104.3 105.9 105.8

0.988 0.978 0.979 0.985 0.982

2.842 2.863 2.880 2.705 2.893

0.529 0.518 0.514 0.573 0.542

6.32 5.24 5.44 4.60 5.02

0.959 0.957

104.2 104.5

0.964

2.912 2.976

0.493

5.02 5.44 ± 0.7

Best ab initio Experiment

The results labelled CRY-x are obtained with the CRYSTAL code, using as basis set 6-31G(d) (x = a), 6-311G(p,d) (x = b), 6-311+G(2d,2p) (x = c); those labelled VASP-a and VASP-b with the VASP code, using ultrasoft pseudo-potentials and projected augmented waves, respectively. In the dimer, dOH and dO–H are the distances of the H atom along the hydrogen bond from the two oxygen atoms, dOO the distance between the oxygens; Ebind are ˚ , angles in °, energies in kcal/mol. For other details, see text. Data in the the binding energies with respect to the isolated molecules. Lengths are in A last two rows are from [32].

114

S. Casassa et al. / Chemical Physics Letters 409 (2005) 110–117

using pseudo-potentials and a plane wave basis set. The same PW91 functional [24] was used as in the CRYSTAL calculations. Core electrons are described by ultra-soft pseudo-potentials [29] and valence electrons by plane wave basis set with a cutoff of 515 eV. This cutoff seems sufficient since the calculated energy of the monomer changes by only 0.2 kcal/mole when passing from 396 to 515 eV. The choice of the pseudo-potentials might play a more important role, and the effect of using projected augmented waves [30] instead has been tested (see below, Section 2.3). Optimizations are carried out with the conjugate gradient algorithm. 2.3. Check of the computational set-up: the water molecule and the water dimer In order to check the adequacy of the present computational set-up, we have considered the performance of the various QM approaches here employed when applied to two simple molecular systems, namely the isolated water molecule and the water dimer. The data for the molecule are anyhow necessary for obtaining ice binding energies. For the VASP calculations, a replicated-molecule model was employed, using a cubic unit ˚ for the molecule and the dimer, cell of side 10 and 15 A respectively. Complete geometry optimization was performed in each case. The results are reported in Table 2. The influence of the GTO basis set on the equilibrium geometry of the water molecule (CRY-x data) is small but not negligible, and for each Hamiltonian there is a certain spread of calculated values. As compared to the experiment, it is seen that HF gives a shorter bond length and a wider angle, while the reverse is observed with DFT-PW91; as expected, the hybrid B3LYP results are intermediate. The VASP result is in good agreement with the corresponding CRYSTAL value as concerns the bond length, but the angle is slightly larger. The case of the dimer is particularly interesting because it is the prototype model of the hydrogen bond which plays an essential role in the stability of ice structures; very accurate QM studies of its structure and properties are available [31,32]. The important data here are the length of the hydroxyl bond of the H atom involved in the hydrogen bond, its ratio with the hydrogen bond,

the distance between the oxygen atoms, which is experimentally accessible, and the energy of formation of the dimer with respect to the isolated molecules. The influence of the quality of the GTO basis set is apparent especially as concerns energy data. For a given Hamiltonian, Ebind is consistently too high with the poorest basis set. The values provided by CRYSTAL with the PW91 Hamiltonian and the two best basis sets are close to the VASP result; the use of projected augmented waves instead of ultra-soft pseudopotentials in VASP would improve the agreement. Geometries are less basis set dependent, while they strongly depend on the Hamiltonian. In particular, HF gives reasonable (but slightly too high) values of the inter-oxygen distance, and too low for the dOH/dO–H ratio; the opposite is true for PW91, while B3LYP gives quite good results in all respects.

3. Results and discussion Table 3 reports energy data for the optimized ice structures, while Tables 4 and 5 give indications about equilibrium geometries, and their departure from ideality. Let us first comment on energies. The order of stability is seen to depend to some extent on the QM approach, although the prevailing indication is that polar structures are more stable than apolar ones, and cubic more stable than hexagonal, contrary to intuition. However, all binding energies obtained with the same technique are coinciding well within 1 kcal/ mol, and their order may be partly related to the degree of optimization reached. The ÔidealÕ values resulting from the VASP-PW91 calculations are even closer among themselves (0.03 kcal/mol); in this respect, it is interesting to note that releasing the condition of ideality is more important for cubic than for hexagonal structures (see below). The difference between the various Hamiltonians is more marked. Pure DFT approaches tend to overestimate, HF to slightly underestimate the experimental value of the binding energy of ordinary ice, which is approximately 14.1 kcal/mol [33], after correction for finite temperature and zero point energy. The overestimation of the hydrogen bond energy in water is a known feature of most

Table 3 Calculated binding energies per H2O (kcal/mol) for the different ice structures and QM approaches referred to the fully optimized water molecule Ice structure

Code and QM Hamiltonian CRYSTAL

PC AC PH AH

VASP

HF

B3LYP

PW91

PW91

10.07 10.39 10.49 10.21

14.49 13.86 14.06 13.88

17.14 16.76 16.93 16.78

17.64 17.50 16.80 16.86

(16.71) (16.70) (16.70) (16.73)

Basis set used for CRYSTAL is 6-311G(p,d). All CRYSTAL data are corrected for BSSE. The last column reports in parenthesis the VASP results for the ideal structures.

S. Casassa et al. / Chemical Physics Letters 409 (2005) 110–117

115

Table 4 Optimized geometrical parameters Ice

Parameter

Code and QM Hamiltonian

Experiment VASP

CRYSTAL HF

B3LYP

PW91

PW91

PC

Volume at ct ˆH HO dOH dO–H

36.71 4.57 7.03 108.0 0.952 1.931

31.61 4.42 6.47 106.5 0.986 1.753

29.41 4.34 6.24 103.9 1.005 1.668

28.76 (29.40) 4.38 (4.36) 6.01 (6.17) 107.5 1.011 (1.008) 1.644 (1.665)

AC

Volume at ct ˆH HO dOH dO–H

36.57 4.79 6.39 106.5 0.952 1.928

31.36 4.50 6.20 109.3 0.987 1.745

29.16 4.38 6.09 107.6 1.007 1.660

28.80 (29.60) 4.27 (4.37) 6.31 (6.19) 107.9 1.009 (1.007) 1.672 (1.672)

PH

Volume a0 b0 c0 ˆH HO dOH dO–H

36.85 4.57 8.40 7.69 108.2 0.952 1.931

31.36 4.42 7.85 7.23 107.4 0.987 1.745

29.33 4.35 7.63 7.07 108.4 1.001 1.665

29.38 (29.37) 4.37 (4.36) 7.56 (7.56) 7.12 (7.12) 109.5 1.008 (1.008) 1.662 (1.664)

32.16 4.5 7.8 7.33 107 0.98 1.800

AH

Volume a0 b0 c0 ˆH HO dOH dO–H

36.71 7.66 8.00 4.79 106.6 0.962 1.854

31.34 7.27 7.66 4.50 107.0 1.000 1.688

29.18 7.08 7.51 4.39 107.6 1.021 1.618

29.42 (29.54) 7.15 (7.14) 7.54 (7.57) 4.37 (4.37) 107.5 1.008 (1.008) 1.666 (1.665)

32.52 7.36 7.36 4.52 107.0 0.985 1.779

˚ 3 per H2O; distances in A ˚ , angles in °. Values for the ÔidealÕ structures obtained with the VASP code are in parenthesis. Experimental Volumes are in A data are taken from [2], Matsuo et al. [Nature (London) 299 (1986) 810] and Leadbetter et al. [J. Chem. Phys. 82 (1985) 424] and concern ice XI and ordinary ice; in the latter case, the lattice parameters of the hcp cell have been ÔtranslatedÕ into those of the corresponding orthorhombic cell.

Table 5 Deformation relative to the parameters for the ideal structures as resulting from VASP-PW91 calculations Ice PC AC PH AH

Parameter ratio pffiffiffi ct =ðat pffiffi2ffiÞ t t c =ða pffiffiffi 2Þ 0 pffiffiffi ðc0 3Þ=ða pffiffiffi 8Þ ðb0 Þ=ða0 3Þ pffiffiffi pffiffiffi 0 8Þ ðc0 3Þ=ða p ffiffi ffi 0 ðb Þ=ða0 3Þ

Calculated value 0.972 1.045 0.998 1.000 1.002 0.997

For ice XI, the two experimental values, to be compared to the calculated ones for PH, are 0.997 and 1.001 Matsuo et al. [Nature (London) 299 (1986) 810] and Leadbetter et al. [J. Chem. Phys. 82 (1985) 424].

DFT techniques, while the standard B3LYP hybrid method provides quite accurate estimates, perhaps due to fortuitous cancellation of errors. This result confirms the indications coming from the study of the water dimer (see Section 2.3). Again, it is found that with the present computational set-up the PW91 binding energies obtained with VASP and CRYSTAL are very close to each other.

As concerns equilibrium geometries, the following can be noted. For all structures, the volume per unit cell is too large with HF, approximately correct with B3LYP, too small with DFT; this takes place with an almost uniform expansion or contraction of the lattice parameters. The essential reason is again to be found in the different description of the hydrogen bond, which is, respectively, too long and too short according to the two kinds of one electron Hamiltonians, while the length of the hydroxyl bond is less dependent on the Hamiltonian adopted. This feature is made more evident when comparing the determinations of the dOH/dO–H ratio with the experimental value, which is 0.544 both in ordinary ice and in ice XI. The average value and the standard deviation of that ratio for the four structures is: 0.493 ± 0.001 (HF); 0.565 ± 0.004 (B3LYP); 0.606 ± 0.008 (PW91); 0.609 ± 0.006 (PW91-VASP). Note incidentally that the length of the two kinds of bond are approximately the same for the ideal as for the optimized structures; however, when the O–H–O is not constrained to remain linear, the dOO distance decreases. The contraction is more pronounced for the cubic structures because of the optimal match between

116

S. Casassa et al. / Chemical Physics Letters 409 (2005) 110–117

the space and the local symmetry: this explains the larger difference in stability between ideal and optimal configurations in these cases. The dependence of angles from the technique adopted is less marked than for lengths. In all cases, in agreement with observation, the angle is wider than in the gas phase molecule. A consequence to be expected from the different description of the hydrogen bond is that proton transfer processes in ice will turn out to be much easier when using DFT rather than HF based models [6–8]. To provide some sort of measure of this difference, we have calculated the energy of the transition state of a hypothetical reaction in cubic ice, where all hydrogen atoms are simultaneously passing from the oxygen to which they are hydroxyl-bonded to the one in front to which they are hydrogen-bonded. The transition state is then a perfect diamond-like structure, with hydrogens midway between oxygens. The optimized O–O distance is practically the same in all cases : 2.39, 2.43, 2.44– ˚ (here and below, the data refer to HF, B3LYP, 2.44 A PW91 in a sequence, a dash separating the CRYSTAL and VASP determination in the last case), corresponding to a strong contraction of the cell volume per water mol˚ 3. ecule in the transition state: 21.0, 22.0, 22.3–22.4 A However, the energy barrier with respect to the equilibrium configuration differs enormously in the three cases: 22.7, 10.7 and 6.1–7.0 kcal/mol, respectively! Consider finally the VASP data reported in Table 5 and concerning departure from ideality of the O sublattice. The maximum departure is observed with the cubic structures, which can be explained by a good match between the local symmetry (C2 axis) and the crystal symmetry (41) in this case: the lattice distortion induces the same distortion of the local tetrahedrons, so allowing the H2O unit to get closer to the molecular geometry. On the contrary, for the hexagonal lattices, in spite of the larger number of independent parameters, the lattice distortion does not improve all the local geometries simultaneously; the resulting geometry, which represents the best compromise, on average, for every local environment, remains closer to the ideal structure. An appreciable contraction in the z direction is observed for the PC structure. The electrostatics is favored by the orientation of the dipoles in the z axis where charges alternate (see Fig. 2): this explains a shortening along the z direction. On the contrary, in the xy plane, the dipoles are parallel and should repel each other. On the whole, the present calculations do not permit a clear-cut discrimination between the different structures on the basis of their stability and of their agreement with the experimental geometric data. Long range forces which favor hexagonal or pseudo-hexagonal structures with respect to cubic-like ones, even in the presence of proton disorder, surely require more sophisticated treatment of the hydrogen bond. This is in fact, as already shown, the weakest aspect of the

QM approaches here adopted. From the point of view of the modeling of ordinary ice, this means that cubiclike or hexagonal-like structures can be chosen as a reference, according to the computational convenience. On the contrary, the use of polar models has to be avoided, as a rule, especially when surface properties are to be simulated with the use of a slab or multi-slab approach [10]. Only apolar structures can generate slabs which are stable, reproduce the bulk properties at their interior when they are sufficiently thick, and present at their surfaces a balanced distribution of dangling hydrogens and dangling oxygens, in order to give the possibility to explore a variety of molecular terminations. Much more critical than structure is the QM approach to be adopted. The present study confirms literature results according to which DFT exchange-correlation potentials in current use underestimate the cost of proton-transfer processes in ice, while the opposite is true for HF. Since at present computer programs for periodic structures are based on one-electron Hamiltonians, this is probably the most serious inconvenience in the use of those tools for the simulation of catalytic processes at ice surfaces. The incorporation of ab initio treatments of electron correlation in periodic programs [34,35] will possibly provide a way out of such difficulties. Acknowledgments This work has been accomplished in the framework of the GDR ÔDynamique Mole´culaire Quantique Applique´e a` la catalyseÕ. Computational facilities provided by IDRIS and CCR are also acknowledged. We are grateful to the referee for his suggestions which have helped us to prepare the present improved version of the Letter, and to Bartolomeo Civalleri for useful discussions and assistance in the computations. Appendix A The complete sets of fractional coordinates for the ÔidealÕ ice structures may be found on the web site: http://www.lct.jussieu.fr/minot/PUBLI/iceweb.pdf. The CRYSTAL03 optimized geometries for the ice structures obtained with the 6311G(p,d) basis set and the different Hamiltonians may be found on the web site: http://www.crystal.unito.it/ looking at ÔApplicationsÕ and ÔSupplementary MaterialÕ. The three basis sets here used, 631G(d), 6311G(p,d) and 6311G+(2p,2d) are also reported there for the sake of reference. References [1] C. Girardet, C. Toubin, Surf. Sci. Rep. 44 (2001) 159. [2] V.F. Petrenko, R.W. Withworth, Physics of Ice, Oxford University Press, 1999.

S. Casassa et al. / Chemical Physics Letters 409 (2005) 110–117 [3] T.K. Hirsch, L. Ojamae, J. Phys. Chem. B 108 (2004) 15856. [4] M. Calatayud, D. Courmier, C. Minot, Chem. Phys. Lett. 369 (2003) 287. [5] J.-L. Kuo, S.J. Singer, Phys. Rev. E 67 (2003) 016114. [6] S. Casassa, C. Pisani, J. Chem. Phys. 116 (2002) 9864. [7] S.C. Xu, J. Chem. Phys. 1R1106 (1999) 2242. [8] Z.F. Liu, C.K. Siu, J.S. Tse, Chem. Phys. Lett. 309 (1999) 335. [9] C. Lee, H. Chen, G. Fitzgerald, J. Chem. Phys. 102 (1995) 1266. [10] G. Bussolin, S. Casassa, C. Pisani, P. Ugliengo, J. Chem. Phys. 108 (1998) 9516. [11] Y. Tajima, T. Matsuo, H. Suga, J. Phys. Chem. Solids 45 (1984) 1135. [12] M.D. Morse, S.A. Rice, J. Chem. Phys. 76 (1982) 650. [13] J.D. Bernal, R.H. Fowler, J. Chem. Phys. 1 (1933) 515. [14] S. Casassa, P. Ugliengo, C. Pisani, J. Chem. Phys. 106 (1997) 8030. [15] P.V. Hobbs, Ice Physics, Clarendon Press, New York, 1974. [16] R. Howe, R.B. Whitworth, J. Chem. Phys. 90 (1989) 4450. [17] I. Minawaga, J. Phys. Soc. Jpn 59 (1990) 1676. [18] S.M. Jackson, V.M. Nield, R.W. Whitworth, M. Oguro, C.C. Wilson, J. Phys. Chem. 101 (1997) 6142. [19] C.M.B. Line, R.W. Whitworth, J. Chem. Phys. 104 (1996) 10008. [20] B.J. Yoon, K. Morokuma, E.R. Davidson, J. Chem. Phys. 83 (1985) 1223. [21] G.T. Barkema, J. de Boer, J. Chem. Phys. 99 (1993) 2059.

117

[22] V.R. Saunders, R. Dovesi, C. Roetti, R. Orlando, C.M. ZicovichWilson, N.M. Harrison, K. Doll, B. Civalleri, I.J. Bush, P. DÕArco, M. LLunell, CRYSTAL 2003 User Manual, Daresbury, Torino, Alessandria, Cuernavaca, London, Braunschweig, Paris, Barcelona, 2003. [23] G. Kresse, J. Hafner, Phys. Rev. B 49 (1994) 14251. [24] J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh, C. Fiolhais, Phys. Rev. B 46 (1992) 6671. [25] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [26] R. Krishnan, J.S. Binkley, R. Seeger, J.A. Pople, J. Chem. Phys. 72 (1980) 650. [27] K. Doll, N.M. Harrison, V.R. Saunders, Int. J. Quantum Chem. 82 (2001) 1. [28] G. Lendvay, I. Mayer, Chem. Phys. Lett. 297 (1998) 365. [29] G. Kresse, J. Hafner, J. Phys. Condens. Matter 6 (1994) 8245. [30] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758. [31] W. Klopper, J.G.C.M. van Duijneveldt-van de Rijdt, F.B. van Duijneveldt, Phys. Chem. Chem. Phys. 2 (2000) 2227. [32] X. Xu, W.A. Goddard III, J. Phys. Chem. A 108 (2004) 2305. [33] E. Whalley, J. Chem. Phys. 81 (1984) 4087. [34] P.Y. Ayala, G.E. Scuseria, J. Chem. Phys. 110 (1999) 3660. [35] C. Pisani, M. Busso, G. Capecchi, S. Casassa, R. Dovesi, L. Maschio, C. Zicovich, M. Schu¨tz, J. Chem. Phys. 122 (2005) 094113.