- Email: [email protected]
Proton ordering energetics in ice phases
Chemical Physics Letters 425 (2006) 246–250 www.elsevier.com/locate/cplett
Proton ordering energetics in ice phases Gareth A. Tribello, Ben Slater
*
Davy Faraday Research Laboratory, The Royal Institution, 21 Albemarle Street, London W1X 4BS, UK Received 24 March 2006; in final form 12 April 2006 Available online 10 May 2006
Abstract Results from first-principles calculations on the subtle energetics of proton ordering in ice phases are shown only to depend on the electrostatic components of the total energy. Proton ordered ice phases can therefore be predicted using electronic structure methods or a tailored potential model. However, analysis of the electron density reveals that high order multipole components, up to hexadecapole, are needed to adequately capture total energy differences between proton ordered and disordered phases. This suggests that current potential models may be unable to reproduce the position of proton ordered ice phases in the phase diagram without extensions to describe high order electrostatics. 2006 Elsevier B.V. All rights reserved.
Solid water ice polymorphs comprise a rich subset of the water phase diagram, which may evolve further as extrema of pressures and temperatures are explored. Currently there are 14 known phases of ice, all of which exhibit tetrahedrally coordinated oxygen atoms. The structural diversity arises because of the multifarious ways that these tetrahedral units can be arranged. The hydrogen atom positions must satisfy the Bernal–Fowler rules [1] that (1) there are two hydrogen atoms adjacent to each oxygen and (2) there is only one hydrogen per hydrogen bond. Pauling [2] showed that for any ice phase there are (3/2)N ways that the hydrogen atoms can be arranged and suggested that all of these potential structures will be degenerate, thus implying all ice structures are proton disordered. However, although proton disordered ice phases occupy the majority of the ice phase diagram, at low temperatures many ice phases undergo a transition to a proton ordered form [3]. In the new ordered phase the underlying oxygen sub-lattice is unchanged but the ordering of the protons gives rise to an increase in symmetry and a change in physicochemical properties. The ordering transition is exemplified by ice Ih (hexagonal ice, the most thermodynamically stable *
Corresponding author. Fax: +44 20 7629 3569. E-mail addresses: [email protected] (G.A. Tribello), [email protected] (B. Slater). 0009-2614/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.04.111
phase of ice at 273 K and ambient pressure), which transforms to ice XI at 72 K and ice VII (a high pressure dense phase of ice) which changes to ice VIII. The structures of ice XI and VIII are shown in Figs. 1a and 1b. There are a number of proton disordered ice phases where the low temperature proton ordered form is as yet unidentified. A longer term aim of this work is to develop a effective potential model that could predict the relative energy of proton networks in ice phases and hence predict the structure and stability of proton ordered and disordered ice. Although a number of interatomic potential models are able to predict the complex features of the water phase diagram [4] with remarkable accuracy, the relative stability of proton disordered and ordered ice phases are usually predicted incorrectly. The absence of a truly transferable water/ice potential that can be used in atomistic simulations remains a considerable challenge to theory. The problem of understanding proton ordering in ice phases is a longstanding one within the literature. Bjerrum [5] was one of the first to approach this phenomenon from a theoretical perspective. He proposed that in the ordered form of ice Ih, the water molecules would be arranged in such a way that suprafacial repulsions between hydrogen atoms on neighbouring water molecules are minimised. If this assertion were true, the ordered form of ice Ih would be anti-ferroelectric, whereas in reality the structure
G.A. Tribello, B. Slater / Chemical Physics Letters 425 (2006) 246–250
Fig. 1a. The structure of ice XI viewed along [1 0 0]. The red spheres indicate oxygen atoms whilst hydrogen is shown in white. The hexagonal symmetry of the oxygen sub-lattice is evident. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 1b. The ice VIII structure viewed along [0 1 0]. The oxygen atoms are shown in red and yellow whilst hydrogen is shown in white. The different coloured oxygen atoms form distinct sub-lattices corresponding to two interpenetrating diamond lattices. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
247
of ice XI is ferroelectric [3]. Furthermore, in the proton ordered ice VIII structure, suprafacial repulsion between hydrogen atoms is not at a minimum. Surprisingly, the vast majority of empirical potential functions predict that the lowest energy structure of ice Ih will be the anti-ferroelectric structure Bjerrum predicted. It has been found by us that even the recently published, extremely sophisticated AMOEBA potential, which features distributed multipoles and Thole type polarisability still predicts the antiferroelectric structure to be more stable than the ferroelectric ice XI phase, irrespective of cutoff and other simulation parameters. However, somewhat reassuringly both Hartree–Fock [6] and periodic density functional theory (DFT) [7] studies correctly predict that the lowest energy hydrogen bonding networks possible for ice Ih and ice VII are the proton ordered XI and VIII structures, respectively. To our knowledge no physical explanation has been given as to why the true structures are lower in energy that the structures Bjerrum and the empirical potentials predict, and clearly this is essential if one wishes to produce and parameterise more accurate potential models that can explain the driving force for proton ordering phase transitions. Recently, Singer et al. [7,8] have developed a technique that uses graph invariants to enumerate symmetry distinct structures of bulk ice and water clusters, and they have shown that the energy of the system (calculated using DFT) can be fitted to a linear combination of these graph invariants. One important conclusion from this work is that the energetics of ice Ih and ice VII can be described by invariants that define the relative orientations of two water molecules connected through only one hydrogen bond. Hirsch and Ojama¨e [9] performed DFT and interatomic potential calculations for 16 symmetrically distinct ice Ih phases and found there was an anti-correlation between the energies calculated by these two methods i.e. low energy proton arrangements are predicted to be high energy using most empirical potentials and vice versa. This observation suggests that a short ranged interaction between neighbouring water molecules is inadequately described by potential models, prompting us to investigate the physical origin of proton ordering energies. In this work, the calculations of Hirsch and Ojama¨e [9] were re-performed on 15 symmetry distinct ice Ih structures using a plane wave approach (CASTEP [10]), a 6 · 3 · 3 Monkhorst-Pack grid of k-points and a plane wave cutoff of 550 eV. In addition, a range of functionals; Perdew-Wang 91 (PW91), Revised Perdew-Burke-Ernzerhof (RPBE) and local density approximation (LDA) were used to explore the sensitivity of relative energy to the functional recipe. All positions and cell parameters are allowed to vary until the energy difference between geometry optimisation steps is less than 5 · 10 6 eV per ion. Fig. 2 gives the total energies, calculated using the PW91, RPBE and LDA functionals, for the 15 symmetry distinct ice Ih configurations relative to the ferroelectric ice XI structure (denoted 0).
248
G.A. Tribello, B. Slater / Chemical Physics Letters 425 (2006) 246–250
Although the absolute cohesive energy might be expected to be very different for each functional recipe, it is clear from Fig. 2 that the three different functionals provide very similar estimates of the energy differences between structures, as has been noted previously [7,9]. Surprisingly, even LDA, which describes the water molecule geometry relatively poorly, gives a similar answer to the two Generalised Gradient Approximation functionals (RPBE, PW91). In the PW91, RPBE and LDA schemes the exchange and correlation terms are treated inequivalently and for many materials, the choice of exchange-correlation recipe strongly affects both the structural and electronic properties. However, in this instance the energy differences are essentially identical within the error limit of these calculations. This suggests the exchange-correlation recipe is not responsible for the observed energy differences and that any tendency of the functional to under- or overbind is cancelled out. By inspection of the DFT total energy expression one infers that the differences in total energy are expected to be due to electrostatic interactions. To substantiate whether the major component of the energy differences between these phases is purely electrostatic in origin, a multipolar analysis was performed. This work exploits the approach developed by Aguado et al. [11] that has been added to the CASTEP code. In essence, the Kohn–Sham orbitals are localised by means of a unitary Wannier transformation [12–14] to form a set of maximally localised Wannier functions (MLWFs). A cutoff radius about every molecule in the system is then defined; any MLWF that lies within the cutoff radius of a particular molecule is then assigned as one of its molecular orbitals. The molecular charge density distribution is then evaluated by summing over these molecular orbitals. Aguado et al. [11] have used this approach for ionic systems and used the resulting charge density to calculate the electrostatic dipoles and quadrupoles at the ion nuclei. In ionic systems,
the dipoles and quadrupoles obtained are those formed by the polarisation of the otherwise ‘spherical’ electron cloud of the ion by the surrounding environment. In this work, the multipoles are expanded about the centre of mass of the molecule and thus incorporate both the static multipoles of the molecule and the induced multipoles. It was decided to use this approach rather than a distributed multipole approach (DMA) because the electron density near the hydrogen atoms will be very small and thus multipoles at the hydrogen atoms will be negligible to the first approximation. In addition, the plane wave method used here prevents us from straightforwardly separating out the atomic contributions to electron density distribution between oxygen and hydrogen atoms. A more sophisticated approach would be to locate multipoles away from the hydrogen nucleus to describe the highly aspherical electron density distribution but the model then becomes more arbitrary and markedly more difficult to implement within, for example, molecular dynamics codes. The interaction energy between the multipoles obtained was then evaluated up to terms in 1/r6 in the electrostatic expansion of the potential energy [15]. Using the 16 symmetry inequivalent structures reported by Hirsch and Ojama¨e for ice Ih and four topologies of ice VII, the multipolar contributions to the total energy differences have been extracted and are displayed in Fig. 3a,b. The first point of interest is that the 1/r3 component of the electrostatic energy (which is that due to dipole–dipole interactions) is strongly anti-correlated with the true total energy differences. Since many interatomic potential models are fitted to reproduce the gas phase dipole, this appears to provide an explanation for why most potential models predict ice XI to be higher in energy than Bjerrum’s phase. The second point of note is that there is excellent agreement between the purely electrostatic description of energy differences and the full DFT energy differences for ice Ih, but
Fig. 2. (a and b) Total energies, relative to the energy of the ice XI structure, for 15 symmetry inequivalent 8 molecule unit cells which represent hypothetical ice Ih structures. The graphs show the correlation between the energies calculated using the PW91 functional and the LDA functional, and the RPBE functional indicating that energy differences between different proton arrangements are essentially independent of the DFT recipe.
G.A. Tribello, B. Slater / Chemical Physics Letters 425 (2006) 246–250
249
Fig. 3. (a) Total energies relative to that of ice XI for the 15 symmetry distinct 8 molecule unit cell configurations of ice Ih. Configuration 0 is the structure of ice XI and configuration 1 is Bjerrum’s predicted ice XI structure. (b) Total energies of 4 symmetry distinct ice VII relative to the energy of ice VIII (configuration 0). In both figures, the dotted line is the PW91 result, filled black circles the electrostatic interaction energy up to terms in 1/r3, empty squares the electrostatic energy up terms in 1/r4, filled black diamonds the electrostatic energy up to terms in 1/r5 and empty triangles denotes the electrostatic energy up to 1/r6. The multipolar energies have been uniformly scaled for ice Ih/XI and VII/VIII to be compatible with the DFT energies, which is necessary because of the absence of a dielectric medium in the electrostatic multipolar calculations.
only when terms up to 1/r6 in the expansion are considered. For ice VII convergence is achieved when terms up 1/r5 are considered. Analyses of the distance dependence of the multipolar components between ice Ih configurations 0 and 1 shows that the dipole–dipole interaction converges very slowly, in agreement with the finding of Buch et al. [16] and Rick [17] who analysed this property for pair potential models. At short range it stabilises configuration 1, whilst at longer range it stabilises configuration 0. As a consequence, the dipole–dipole contribution to the energy difference between structures is actually rather small, so that higher order multipolar terms become important. The rapid decay of 1/r4, 1/ r5 and 1/r6 terms ensures that higher order terms converge by the second coordination sphere. The orientations of molecules in the second coordination are a geometric consequence of the topology of the first coordination sphere, which explains why the invariant approach of Singer et al., which fits the energy to very local geometric features, is able to correctly describe variation in energy as the protons are permuted. The driving force for proton ordering energetics in ice phases can thus be understood as a purely electrostatic phenomenon. Note that this result is distinct from observations that hydrogen bonding is primarily due to electrostatic terms. What has been observed here is that structures with identical oxygen sub-lattices but with different hydrogen bonding networks have distinct energies and that the energy differences are only due to variance in the electrostatic interactions. Furthermore, it is unsurprising that potentials, which are generally fitted to the gas phase electrostatics of the water molecule, fail to describe the energetics adequately. The electronic cloud around the water molecule is distorted in ice and it will be necessary to account for this distortion when fitting potentials to reproduce proton ordering. Since proton ordering is purely driven by the electrostatic field, it follows that the transition temperatures will
be strongly influenced by local electrical fields within the crystallite due to, for example, lattice impurities, which explains why moderate variation in experimentally determined transition temperatures has been observed. The surprisingly high order multipoles required to describe the electrostatic energy highlights a serious deficiency in most literature empirical which are typically fitted to the gas phase dipole and polarisability. Using distributed multipoles it may be possible to obtain a reasonable description of the charge density using fewer terms. However, the magnitude of the multipoles in both the DMA and the molecule centred multipole scheme will be dominated by the oxygen atom and hence from the observations presented here, it is probable that a similar order of multipole will be necessary within a DMA. These results indicate that a successful ice potential must describe electronic deformation which results in the formation of induced quadrupoles, octupoles etc., which would be both computationally expensive and difficult to implement. An alternative would be to derive a potential, which contains a term or terms that mimic the effect of these high order electrostatics. For transferability, we note that a potential would have to account for complex deformations of the electron density and hence a polarisable model that builds upon the static multipole is almost certainly necessary. These findings support and build upon those of Xantheas et al. [18,19] who showed that terms up to the hexadecapoles are necessary to describe the electric fields in bulk ice and water clusters. Given that the vast majority of forcefields are fitted only to reproduce dipolar properties it is unsurprising that they predict the incorrect proton arrangement and relative energy for ice XI (see Buch et al. [16] and Rick [17] for related discussion). The significance of the proton ordering transitions and the ability of potential models to predict the correct structure/energy relationships pervades beyond ice phases. At interfaces between water and sorbate molecules, water is known to
250
G.A. Tribello, B. Slater / Chemical Physics Letters 425 (2006) 246–250
order [20,21] and it is therefore conceivable that literature potentials bias proton ordering networks towards high energy anti-ferroelectric arrangements rather than low energy ferroelectric arrangements. As a consequence, delicate processes such as free-energy barriers between protein conformers and barriers to adsorption on ordered or disordered solid materials may be affected. In addition, water clustering in confined hosts such as zeolites and other nanoporous hosts may be subtly perturbed by deficiencies within existing potentials. Development of fully transferable water potentials capable of describing the ice phase diagram should therefore incorporate or mimic multipolar terms in order to reproduce the subtle proton ordering energetics. Finally, since relatively standard periodic electronic structure codes have the intrinsic accuracy to resolve between proton ordered and disordered states, it follows that the total energy from, for example, DFT codes can be used as a predictive tool to determine the proton ordered states of known proton disordered ice phases. Acknowledgments Thanks to Prof. Paul Madden, Dr. Andre´s Aguado and Lindsay Foy for help and discussions on the Wannier function calculations. Thanks also to EPSRC for an allocation of computing time on the HPC(x) supercomputing facility via the Materials Chemistry consortium.
References [1] J.D. Bernal, R.H. Fowler, J. Chem. Phys. 1 (1933) 515. [2] L. Pauling, J. Am. Chem. Soc. 57 (1935) 2680. [3] V.F. Petrenko, R.W. Whitworth, Physics of Ice, Oxford University Press, 2002. [4] E. Sanz, C. Vega, J.L. Abascal, L.G. MacDowell, J. Chem. Phys. 121 (2004) 1165. [5] N. Bjerrum, Science 115 (1952) 386. [6] S. Casassa, M. Calatayud, K. Doll, C. Minot, C. Pisani, Chem. Phys. Lett. 409 (2005) 110. [7] S.J. Singer, J.-L. Kuo, T.K. Hirsch, C. Knight, L. Ojama¨e, M.L. Klein, Phys. Rev. Lett. 94 (2005) 135701. [8] J.-L. Kuo, S.J. Singer, Phys. Rev. E 67 (2003) 016114. [9] T.K. Hirsch, L. Ojama¨e, J. Phys. Chem. B 108 (2004) 15856. [10] M.D. Segall, P.J.D. Lindan, M.J. Probert, C.J. Pickard, P.J. Hasnip, S.J. Clark, M.C. Payne, J. Phys.: Cond. Matt. 14 (2002) 2712. [11] A. Aguado, L. Bernasconi, S. Jahn, P.A. Madden, Faraday Discuss. 124 (2003) 171. [12] N.W. Ashcroft, N.D. Mermin, Solid State Phys. (1976). [13] N. Mazari, D. Vanderbilt, Phys. Rev. B 56 (1997) 12847. [14] E.I. Blout, Solid State Phys. 13 (1962) 305. [15] A.J. Stone, Theory of Intermolecular Forces, Oxford University Press, 2000. [16] V. Buch, P. Sandler, J. Sadlej, Phys. Chem. B 102 (1998) 8641. [17] S.W. Rick, J. Chem. Phys. 122 (2005) 094504. [18] E.R. Batista, S.S. Xantheas, H. Jonson, J. Chem. Phys. 109 (1998) 4546. [19] E.R. Batista, S.S. Xantheas, H. Jonson, J. Chem. Phys. 112 (2000) 3285. [20] S. Keresit, S.C. Parker, J. Am. Chem. Soc. 126 (2004) 10152. [21] C.Y. Ruan, V.A. Lobastov, F. Vigliotti, S.Y. Chen, A.H. Zewail, Science 304 (2004) 80.
Proton ordering energetics in ice phases Gareth A. Tribello, Ben Slater
*
Davy Faraday Research Laboratory, The Royal Institution, 21 Albemarle Street, London W1X 4BS, UK Received 24 March 2006; in final form 12 April 2006 Available online 10 May 2006
Abstract Results from first-principles calculations on the subtle energetics of proton ordering in ice phases are shown only to depend on the electrostatic components of the total energy. Proton ordered ice phases can therefore be predicted using electronic structure methods or a tailored potential model. However, analysis of the electron density reveals that high order multipole components, up to hexadecapole, are needed to adequately capture total energy differences between proton ordered and disordered phases. This suggests that current potential models may be unable to reproduce the position of proton ordered ice phases in the phase diagram without extensions to describe high order electrostatics. 2006 Elsevier B.V. All rights reserved.
Solid water ice polymorphs comprise a rich subset of the water phase diagram, which may evolve further as extrema of pressures and temperatures are explored. Currently there are 14 known phases of ice, all of which exhibit tetrahedrally coordinated oxygen atoms. The structural diversity arises because of the multifarious ways that these tetrahedral units can be arranged. The hydrogen atom positions must satisfy the Bernal–Fowler rules [1] that (1) there are two hydrogen atoms adjacent to each oxygen and (2) there is only one hydrogen per hydrogen bond. Pauling [2] showed that for any ice phase there are (3/2)N ways that the hydrogen atoms can be arranged and suggested that all of these potential structures will be degenerate, thus implying all ice structures are proton disordered. However, although proton disordered ice phases occupy the majority of the ice phase diagram, at low temperatures many ice phases undergo a transition to a proton ordered form [3]. In the new ordered phase the underlying oxygen sub-lattice is unchanged but the ordering of the protons gives rise to an increase in symmetry and a change in physicochemical properties. The ordering transition is exemplified by ice Ih (hexagonal ice, the most thermodynamically stable *
Corresponding author. Fax: +44 20 7629 3569. E-mail addresses: [email protected] (G.A. Tribello), [email protected] (B. Slater). 0009-2614/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.04.111
phase of ice at 273 K and ambient pressure), which transforms to ice XI at 72 K and ice VII (a high pressure dense phase of ice) which changes to ice VIII. The structures of ice XI and VIII are shown in Figs. 1a and 1b. There are a number of proton disordered ice phases where the low temperature proton ordered form is as yet unidentified. A longer term aim of this work is to develop a effective potential model that could predict the relative energy of proton networks in ice phases and hence predict the structure and stability of proton ordered and disordered ice. Although a number of interatomic potential models are able to predict the complex features of the water phase diagram [4] with remarkable accuracy, the relative stability of proton disordered and ordered ice phases are usually predicted incorrectly. The absence of a truly transferable water/ice potential that can be used in atomistic simulations remains a considerable challenge to theory. The problem of understanding proton ordering in ice phases is a longstanding one within the literature. Bjerrum [5] was one of the first to approach this phenomenon from a theoretical perspective. He proposed that in the ordered form of ice Ih, the water molecules would be arranged in such a way that suprafacial repulsions between hydrogen atoms on neighbouring water molecules are minimised. If this assertion were true, the ordered form of ice Ih would be anti-ferroelectric, whereas in reality the structure
G.A. Tribello, B. Slater / Chemical Physics Letters 425 (2006) 246–250
Fig. 1a. The structure of ice XI viewed along [1 0 0]. The red spheres indicate oxygen atoms whilst hydrogen is shown in white. The hexagonal symmetry of the oxygen sub-lattice is evident. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 1b. The ice VIII structure viewed along [0 1 0]. The oxygen atoms are shown in red and yellow whilst hydrogen is shown in white. The different coloured oxygen atoms form distinct sub-lattices corresponding to two interpenetrating diamond lattices. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
247
of ice XI is ferroelectric [3]. Furthermore, in the proton ordered ice VIII structure, suprafacial repulsion between hydrogen atoms is not at a minimum. Surprisingly, the vast majority of empirical potential functions predict that the lowest energy structure of ice Ih will be the anti-ferroelectric structure Bjerrum predicted. It has been found by us that even the recently published, extremely sophisticated AMOEBA potential, which features distributed multipoles and Thole type polarisability still predicts the antiferroelectric structure to be more stable than the ferroelectric ice XI phase, irrespective of cutoff and other simulation parameters. However, somewhat reassuringly both Hartree–Fock [6] and periodic density functional theory (DFT) [7] studies correctly predict that the lowest energy hydrogen bonding networks possible for ice Ih and ice VII are the proton ordered XI and VIII structures, respectively. To our knowledge no physical explanation has been given as to why the true structures are lower in energy that the structures Bjerrum and the empirical potentials predict, and clearly this is essential if one wishes to produce and parameterise more accurate potential models that can explain the driving force for proton ordering phase transitions. Recently, Singer et al. [7,8] have developed a technique that uses graph invariants to enumerate symmetry distinct structures of bulk ice and water clusters, and they have shown that the energy of the system (calculated using DFT) can be fitted to a linear combination of these graph invariants. One important conclusion from this work is that the energetics of ice Ih and ice VII can be described by invariants that define the relative orientations of two water molecules connected through only one hydrogen bond. Hirsch and Ojama¨e [9] performed DFT and interatomic potential calculations for 16 symmetrically distinct ice Ih phases and found there was an anti-correlation between the energies calculated by these two methods i.e. low energy proton arrangements are predicted to be high energy using most empirical potentials and vice versa. This observation suggests that a short ranged interaction between neighbouring water molecules is inadequately described by potential models, prompting us to investigate the physical origin of proton ordering energies. In this work, the calculations of Hirsch and Ojama¨e [9] were re-performed on 15 symmetry distinct ice Ih structures using a plane wave approach (CASTEP [10]), a 6 · 3 · 3 Monkhorst-Pack grid of k-points and a plane wave cutoff of 550 eV. In addition, a range of functionals; Perdew-Wang 91 (PW91), Revised Perdew-Burke-Ernzerhof (RPBE) and local density approximation (LDA) were used to explore the sensitivity of relative energy to the functional recipe. All positions and cell parameters are allowed to vary until the energy difference between geometry optimisation steps is less than 5 · 10 6 eV per ion. Fig. 2 gives the total energies, calculated using the PW91, RPBE and LDA functionals, for the 15 symmetry distinct ice Ih configurations relative to the ferroelectric ice XI structure (denoted 0).
248
G.A. Tribello, B. Slater / Chemical Physics Letters 425 (2006) 246–250
Although the absolute cohesive energy might be expected to be very different for each functional recipe, it is clear from Fig. 2 that the three different functionals provide very similar estimates of the energy differences between structures, as has been noted previously [7,9]. Surprisingly, even LDA, which describes the water molecule geometry relatively poorly, gives a similar answer to the two Generalised Gradient Approximation functionals (RPBE, PW91). In the PW91, RPBE and LDA schemes the exchange and correlation terms are treated inequivalently and for many materials, the choice of exchange-correlation recipe strongly affects both the structural and electronic properties. However, in this instance the energy differences are essentially identical within the error limit of these calculations. This suggests the exchange-correlation recipe is not responsible for the observed energy differences and that any tendency of the functional to under- or overbind is cancelled out. By inspection of the DFT total energy expression one infers that the differences in total energy are expected to be due to electrostatic interactions. To substantiate whether the major component of the energy differences between these phases is purely electrostatic in origin, a multipolar analysis was performed. This work exploits the approach developed by Aguado et al. [11] that has been added to the CASTEP code. In essence, the Kohn–Sham orbitals are localised by means of a unitary Wannier transformation [12–14] to form a set of maximally localised Wannier functions (MLWFs). A cutoff radius about every molecule in the system is then defined; any MLWF that lies within the cutoff radius of a particular molecule is then assigned as one of its molecular orbitals. The molecular charge density distribution is then evaluated by summing over these molecular orbitals. Aguado et al. [11] have used this approach for ionic systems and used the resulting charge density to calculate the electrostatic dipoles and quadrupoles at the ion nuclei. In ionic systems,
the dipoles and quadrupoles obtained are those formed by the polarisation of the otherwise ‘spherical’ electron cloud of the ion by the surrounding environment. In this work, the multipoles are expanded about the centre of mass of the molecule and thus incorporate both the static multipoles of the molecule and the induced multipoles. It was decided to use this approach rather than a distributed multipole approach (DMA) because the electron density near the hydrogen atoms will be very small and thus multipoles at the hydrogen atoms will be negligible to the first approximation. In addition, the plane wave method used here prevents us from straightforwardly separating out the atomic contributions to electron density distribution between oxygen and hydrogen atoms. A more sophisticated approach would be to locate multipoles away from the hydrogen nucleus to describe the highly aspherical electron density distribution but the model then becomes more arbitrary and markedly more difficult to implement within, for example, molecular dynamics codes. The interaction energy between the multipoles obtained was then evaluated up to terms in 1/r6 in the electrostatic expansion of the potential energy [15]. Using the 16 symmetry inequivalent structures reported by Hirsch and Ojama¨e for ice Ih and four topologies of ice VII, the multipolar contributions to the total energy differences have been extracted and are displayed in Fig. 3a,b. The first point of interest is that the 1/r3 component of the electrostatic energy (which is that due to dipole–dipole interactions) is strongly anti-correlated with the true total energy differences. Since many interatomic potential models are fitted to reproduce the gas phase dipole, this appears to provide an explanation for why most potential models predict ice XI to be higher in energy than Bjerrum’s phase. The second point of note is that there is excellent agreement between the purely electrostatic description of energy differences and the full DFT energy differences for ice Ih, but
Fig. 2. (a and b) Total energies, relative to the energy of the ice XI structure, for 15 symmetry inequivalent 8 molecule unit cells which represent hypothetical ice Ih structures. The graphs show the correlation between the energies calculated using the PW91 functional and the LDA functional, and the RPBE functional indicating that energy differences between different proton arrangements are essentially independent of the DFT recipe.
G.A. Tribello, B. Slater / Chemical Physics Letters 425 (2006) 246–250
249
Fig. 3. (a) Total energies relative to that of ice XI for the 15 symmetry distinct 8 molecule unit cell configurations of ice Ih. Configuration 0 is the structure of ice XI and configuration 1 is Bjerrum’s predicted ice XI structure. (b) Total energies of 4 symmetry distinct ice VII relative to the energy of ice VIII (configuration 0). In both figures, the dotted line is the PW91 result, filled black circles the electrostatic interaction energy up to terms in 1/r3, empty squares the electrostatic energy up terms in 1/r4, filled black diamonds the electrostatic energy up to terms in 1/r5 and empty triangles denotes the electrostatic energy up to 1/r6. The multipolar energies have been uniformly scaled for ice Ih/XI and VII/VIII to be compatible with the DFT energies, which is necessary because of the absence of a dielectric medium in the electrostatic multipolar calculations.
only when terms up to 1/r6 in the expansion are considered. For ice VII convergence is achieved when terms up 1/r5 are considered. Analyses of the distance dependence of the multipolar components between ice Ih configurations 0 and 1 shows that the dipole–dipole interaction converges very slowly, in agreement with the finding of Buch et al. [16] and Rick [17] who analysed this property for pair potential models. At short range it stabilises configuration 1, whilst at longer range it stabilises configuration 0. As a consequence, the dipole–dipole contribution to the energy difference between structures is actually rather small, so that higher order multipolar terms become important. The rapid decay of 1/r4, 1/ r5 and 1/r6 terms ensures that higher order terms converge by the second coordination sphere. The orientations of molecules in the second coordination are a geometric consequence of the topology of the first coordination sphere, which explains why the invariant approach of Singer et al., which fits the energy to very local geometric features, is able to correctly describe variation in energy as the protons are permuted. The driving force for proton ordering energetics in ice phases can thus be understood as a purely electrostatic phenomenon. Note that this result is distinct from observations that hydrogen bonding is primarily due to electrostatic terms. What has been observed here is that structures with identical oxygen sub-lattices but with different hydrogen bonding networks have distinct energies and that the energy differences are only due to variance in the electrostatic interactions. Furthermore, it is unsurprising that potentials, which are generally fitted to the gas phase electrostatics of the water molecule, fail to describe the energetics adequately. The electronic cloud around the water molecule is distorted in ice and it will be necessary to account for this distortion when fitting potentials to reproduce proton ordering. Since proton ordering is purely driven by the electrostatic field, it follows that the transition temperatures will
be strongly influenced by local electrical fields within the crystallite due to, for example, lattice impurities, which explains why moderate variation in experimentally determined transition temperatures has been observed. The surprisingly high order multipoles required to describe the electrostatic energy highlights a serious deficiency in most literature empirical which are typically fitted to the gas phase dipole and polarisability. Using distributed multipoles it may be possible to obtain a reasonable description of the charge density using fewer terms. However, the magnitude of the multipoles in both the DMA and the molecule centred multipole scheme will be dominated by the oxygen atom and hence from the observations presented here, it is probable that a similar order of multipole will be necessary within a DMA. These results indicate that a successful ice potential must describe electronic deformation which results in the formation of induced quadrupoles, octupoles etc., which would be both computationally expensive and difficult to implement. An alternative would be to derive a potential, which contains a term or terms that mimic the effect of these high order electrostatics. For transferability, we note that a potential would have to account for complex deformations of the electron density and hence a polarisable model that builds upon the static multipole is almost certainly necessary. These findings support and build upon those of Xantheas et al. [18,19] who showed that terms up to the hexadecapoles are necessary to describe the electric fields in bulk ice and water clusters. Given that the vast majority of forcefields are fitted only to reproduce dipolar properties it is unsurprising that they predict the incorrect proton arrangement and relative energy for ice XI (see Buch et al. [16] and Rick [17] for related discussion). The significance of the proton ordering transitions and the ability of potential models to predict the correct structure/energy relationships pervades beyond ice phases. At interfaces between water and sorbate molecules, water is known to
250
G.A. Tribello, B. Slater / Chemical Physics Letters 425 (2006) 246–250
order [20,21] and it is therefore conceivable that literature potentials bias proton ordering networks towards high energy anti-ferroelectric arrangements rather than low energy ferroelectric arrangements. As a consequence, delicate processes such as free-energy barriers between protein conformers and barriers to adsorption on ordered or disordered solid materials may be affected. In addition, water clustering in confined hosts such as zeolites and other nanoporous hosts may be subtly perturbed by deficiencies within existing potentials. Development of fully transferable water potentials capable of describing the ice phase diagram should therefore incorporate or mimic multipolar terms in order to reproduce the subtle proton ordering energetics. Finally, since relatively standard periodic electronic structure codes have the intrinsic accuracy to resolve between proton ordered and disordered states, it follows that the total energy from, for example, DFT codes can be used as a predictive tool to determine the proton ordered states of known proton disordered ice phases. Acknowledgments Thanks to Prof. Paul Madden, Dr. Andre´s Aguado and Lindsay Foy for help and discussions on the Wannier function calculations. Thanks also to EPSRC for an allocation of computing time on the HPC(x) supercomputing facility via the Materials Chemistry consortium.
References [1] J.D. Bernal, R.H. Fowler, J. Chem. Phys. 1 (1933) 515. [2] L. Pauling, J. Am. Chem. Soc. 57 (1935) 2680. [3] V.F. Petrenko, R.W. Whitworth, Physics of Ice, Oxford University Press, 2002. [4] E. Sanz, C. Vega, J.L. Abascal, L.G. MacDowell, J. Chem. Phys. 121 (2004) 1165. [5] N. Bjerrum, Science 115 (1952) 386. [6] S. Casassa, M. Calatayud, K. Doll, C. Minot, C. Pisani, Chem. Phys. Lett. 409 (2005) 110. [7] S.J. Singer, J.-L. Kuo, T.K. Hirsch, C. Knight, L. Ojama¨e, M.L. Klein, Phys. Rev. Lett. 94 (2005) 135701. [8] J.-L. Kuo, S.J. Singer, Phys. Rev. E 67 (2003) 016114. [9] T.K. Hirsch, L. Ojama¨e, J. Phys. Chem. B 108 (2004) 15856. [10] M.D. Segall, P.J.D. Lindan, M.J. Probert, C.J. Pickard, P.J. Hasnip, S.J. Clark, M.C. Payne, J. Phys.: Cond. Matt. 14 (2002) 2712. [11] A. Aguado, L. Bernasconi, S. Jahn, P.A. Madden, Faraday Discuss. 124 (2003) 171. [12] N.W. Ashcroft, N.D. Mermin, Solid State Phys. (1976). [13] N. Mazari, D. Vanderbilt, Phys. Rev. B 56 (1997) 12847. [14] E.I. Blout, Solid State Phys. 13 (1962) 305. [15] A.J. Stone, Theory of Intermolecular Forces, Oxford University Press, 2000. [16] V. Buch, P. Sandler, J. Sadlej, Phys. Chem. B 102 (1998) 8641. [17] S.W. Rick, J. Chem. Phys. 122 (2005) 094504. [18] E.R. Batista, S.S. Xantheas, H. Jonson, J. Chem. Phys. 109 (1998) 4546. [19] E.R. Batista, S.S. Xantheas, H. Jonson, J. Chem. Phys. 112 (2000) 3285. [20] S. Keresit, S.C. Parker, J. Am. Chem. Soc. 126 (2004) 10152. [21] C.Y. Ruan, V.A. Lobastov, F. Vigliotti, S.Y. Chen, A.H. Zewail, Science 304 (2004) 80.