Periodic MP2, RPA, and Boundary Condition Assessment of Hydrogen

Letter pubs.acs.org/JPCL

Periodic MP2, RPA, and Boundary Condition Assessment of Hydrogen Ordering in Ice XV Mauro Del Ben,† Joost VandeVondele,*,‡ and Ben Slater*,§ †

Department of Chemistry, University of Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland Department of Materials, ETH Zürich, Wolfgang-Pauli-Strasse 27, CH-8093 Zürich, Switzerland § Department of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom ‡

S Supporting Information *

ABSTRACT: Ice XV is the hydrogen-ordered form of the ice VI phase whose structure was predicted to be Cc and ferroelectric using periodic DFT approaches. However, neutron diffraction and Raman spectroscopy data show the structure to have P1̅ symmetry and to be antiferroelectric. Recent work1 using fragment-based MP2 and CCSD(T) approaches predicts the experimental structure as the ground state. We have reconsidered this problem using fully periodic MP2 and RPA approaches and find that the ferroelectric Cc structure is the lowest energy configuration. However, ubiquitously employed tinfoil boundary conditions stabilize polar structures. We suggest that ferroelectric Cc crystals can grow within a paraelectric ice VI matrix but may become unstable once a fraction of the matrix has become hydrogen-ordered. The reduction in dielectric constant causes P1̅ and other structures with small polarization to become favored, providing a possible resolution between observation and theoretical predictions. SECTION: Molecular Structure, Quantum Chemistry, and General Theory

I

ce exhibits 15 distinct crystalline phases,2 and the chemistry of the most readily accessible phase, ice Ih, continues to yield surprises, such as “cubic ice”, which is more accurately described as stacking disordered ice Isd,3 a newly proposed ice 0 phase4−6 that appears to facilitate the homogeneous nucleation of ice and the recent experimental identification of supercooled water beyond its accepted homogeneous freezing point.7 High-pressure ice phases have also attracted attention, and the most recently identified stable form of ice, ice XV,8 is a case in point. Almost a decade ago, Knight and Singer9 and Kuo and Kuhs10,11 independently predicted that ice XV, a stable hydrogen-disordered ice phase in the ∼1 GPa regime,should have a crystallographically distinct hydrogen-ordered structural form of Cc space group symmetry. Both groups predicted the same crystal structure, obtained using local density functional approaches, often referred to in the literature2,8 as the KSKK structure. When Salzmann et al. isolated and reported the hydrogen ordered form of ice VI in 2009, it was found to have the P1̅ crystal structure in contrast with KSKK prediction. The discrepancy between theory and experiment was surprising in that numerous density functional theory studies at various levels of sophistication had successfully identified the hydrogenordered counterparts of ice I (ice XI12−14), ice II,15 ice III (ice IX15), ice V (ice XII16), and ice VII(ice VIII12,17,18). See also the review by Singer and Knight.19 The mismatch between theory and experiment was highlighted in a 2011 review2 where ice XV was proposed as a solid-state benchmark system for hydrogen bonding. In 2013, Nanda and Beran1 rose to the challenge and by using a hierarchical, fragment-based quantum mechanical/molecular mechanics approach, they found that the © XXXX American Chemical Society

experimentally identified structure was favored over other structures at the MP2 and CCSD(T) levels, apparently resolving the conflict between experiment and theory. Aside from the crystallographic differences and corresponding electrical properties of each structure, another point of interest is that ice XV has the largest variation in hydrogen bonding angles of any crystalline ice phase. The latter property suggests that ice XV is a testing material for theoretical approaches; modeling of hydrogen bonds under pressure suggests that currently no density functional exists that can consistently predict density and the cohesive energies of phases with an accuracy comparable to diffusion Monte Carlo results.20,21 Moreover, in ice XV there is a coupled problem of describing the highly distorted arrangement of hydrogen bonds within each sublattice and the separation between the independent sublattices, which is affected by van der Waals interactions between the sublattices. van der Waals interactions are currently an area of very active development22,23 in the density functional field, and in the absence of functionals that can provide chemical accuracy for solid cohesive energies over a wide range of densities there has been progress in implementing embedding techniques,24,25 periodic HF with CCSD corrections,26 periodic CCSD,27 local periodic MP2 methods,28,29 and random phase approximation (RPA) approaches.30 This motivated us to examine the ice XV Received: September 18, 2014 Accepted: November 11, 2014

4122

dx.doi.org/10.1021/jz501985w | J. Phys. Chem. Lett. 2014, 5, 4122−4128

The Journal of Physical Chemistry Letters

Letter

the experimentally identified structure (which we denote SRMF) and in Figure 1b we show the structure predicted by Knight and Singer and Kuo and Kuhs (KSKK). There is an important difference in the two lattices: the SRMF structure is antiferroelectric and the KSKK structure is ferroelectric. The structure consists of two interpenetrating lattices with the Edingtonite topology;35 in the SRMF structure each sublattice is polar (and maximally so2), but the inversion symmetry of the P1̅ space group means that the dipole moments of each sublattice cancel. In the KSKK structure each sublattice is also polar and the dipole moments of each sublattice only partially cancel, giving rise to an overall net dipole for the unit cell. Using the oxygen net of ice VI, an individual network can hydrogen order in three different ways, which have been labeled A−C, in accord with the original report of ice XV.8 The unit cell contains 10 water molecules, and according to symmetry analysis (using graph invariant software36), there are 18 symmetry inequivalent structures within the space groups 2, 4, 7, and 9 (corresponding to P1̅, P21, Pn, and Cc, respectively).2 Note that there are two nonequivalent structures in the Cc (9) space group, which are labeled 9A1 and 9A2 after ref 1. According to our analysis, 9A2 (see Figure 1b) is isostructural with the KSKK structure, whereas ref 1 reports that the 9A1 (see Figure 1c) structure is the KSKK structure. 9A1 can be transformed to 9A2 by a 180° rotation of one of its sublattices about the c axis. The structural difference between 9A1 and 9A2 has a profound effect on their polarity, as will be seen later. Next, as a consistency check for our computational methodology, we performed a full cell optimization (without symmetry constraints, i.e., P1 symmetry) on all 18 of the symmetry inequivalent structures at the GGA PBE37 level of theory (which yields reliable structures, see, e.g., ref 38) using two different computer programs. In Figure 2, the total energies are reported relative to the average energy computed using CP2K, which uses a dual local orbital and plane-wave-based scheme and CASTEP,39 which uses different pseudopotentials and also a plane-wave basis and is hence basis-set-superposition

problem using fully periodic RPA and MP2 treatments that have been recently implemented31,32 in the CP2K/Quickstep code.33,34 Although several studies of ice phases28,30 have been reported at the periodic MP2 and RPA levels, this is the first application of these methods to ice XV. We begin by very briefly reviewing the structure of ice XV; in Figure 1a, we show

Figure 1. Three of the 18 symmetry distinct structures that are possible within the ice XV unit cell. Hydrogen is shown as white spheres. (a) Experimentally determined SRMF ice XV structure (2C1) exhibiting P1̅ symmetry. (b) KSKK (9A2) structure predicted by firstprinciples approaches exhibiting Cc symmetry. (c) Other Cc symmetry structure (9A1) that has a distinct hydrogen-bonding pattern and sublattice orientation from the KSKK structure. Oxygen atoms are shown in green and red to distinguish the distinct sublattices within the structure. In each case a 2 × 2 × 2 supercell is shown along the [001] axis.

Figure 2. Relative energies, with respect to the average, of the 18 symmetry inequivalent structures possible in the ice XV unit cell reported at the PBE level of theory, consistently optimized with both codes. Energies on the vertical axis are expressed in kilojoules per mole per molecule relative to the average energy of all 18 structures and for the unit cell containing 10 molecules. Each structure type is labeled after the convention of ref 1. 4123



Letter

Figure 3. Relative energies in kilojoules per mole per molecule, with respect to the average, of the 18 symmetry inequivalent structures possible in the ice XV unit cell reported at the PBE, PBE0, and PBE0-D3 level of theory employing MP2-optimized geometries.

Figure 4. Relative energies in kilojoules per mole per molecule, with respect to the average, of the 18 symmetry inequivalent structures possible in the ice XV unit cell reported at the periodic RI-MP2 and RI-RPA level of theory using the cc-QZ basis employing MP2-optimized geometries.

the band gap, which is, for example, important to estimate the static dielectric constant of ice Ih.41 In the latter scheme, D3 indicates the use of Grimme’s van der Waals correction,42 which takes into account interactions that are missing from GGA density functionals in a semiempirical way. The results are shown in Figure 3 and are remarkably consistent with the PBE results, and, in particular, the 9A2 structure is again found to be the global minimum. Although relative energies of all structures are insensitive to the functional scheme employed, there is some variation in the absolute density. PBE gives 1.40 g/cm3, PBE0 1.40 to 1.41 g/cm3, and PBE0-D3 1.47 to 1.48 g/cm3, in accord with similar studies of ice phases20,21 and liquid water43−45 Given the wide variations in intramolecular bond angle, variation in intermolecular hydrogen bonding angle, and nonequivalent van der Waals interaction between the sublattices, the relative energies of the 18 possible structures could be sensitive to the precise formulation of the density

error-free. Details of all simulation settings used in this letter and a more detailed study of basis set convergence are described in the Supporting Information (SI). The largest difference in energy was found to be ∼0.03 kJ/mol and suggests that the two computational procedures have similar accuracy and also validates the use of atom-centered orbitals, even for these small relative energies. In agreement with previous works,9−11 it was found that the global minimum is the KSKK structure (Cc space group), which is structure 9A2 in Figure 2. The experimentally identified structure is 2C1 and is unambiguously metastable with respect to several A, B, and C type structures. To investigate if variations in exchange and correlation functional have a significant impact, we examine what the hybrid PBE0 functional,40 with 25% HF exchange used in place of PBE exchange, as well as what the PBE0-D3 scheme predicts. In the former scheme, the introduction of HF exchange opens 4124



Letter

functional employed. However, we find this is not the case for ice XV. In the Supporting Information (SI) we present more data to show that this holds more broadly. Summarizing the SI, other approaches for incorporating van der Waals interaction, based on a nonlocal kernel and the electron density, such as vdw-DF2,46 optPBE-B88vdw,47 and rVV10,48 also show very similar relative energies and also give 9A2 as the minimum energy structure. Also, changing from the PBE family of functionals to the M06 family49 of functionals, with and without dispersion corrections, and including M06-2X-D3 scheme which uses 54% of HF exchange, leads to very similar results. A similar observation has been made previously for I/XI,50 which has a much smaller variance in the intermolecular hydrogen bonding angles. We also show that the zero-point vibrational energy contributions do not significantly affect the relative energies of the 18 possible structures in the ice XV unit cell. Higher level methods including wave-function correlation offer an independent way of assessing the veracity of the density functional predictions. Nanda and Beran’s work based on MP2 and CCSD(T), using a fragment-based, hybrid manybody interaction, quantum-mechanical/molecular mechanics method, energetically favors the experimentally resolved P1̅ (2C1) structure over the 2A1, 2B1, 9A1 and 9B2 structures. At the MP2 level, the remaining 13 structures were found to be significantly less stable. In this work, we also examine all 18 structures using fully periodic MP2 and RPA approaches. RPA can be viewed as an approximate CCSD method that only retains part of the diagrams.51 RPA has been shown to yield relative energies for ice phases that approach the accuracy of quantum Monte Carlo methods.30 For the MP2 method, forces and stresses have been implemented in CP2K, allowing fully consistent energies and geometries for a high level of theory. Basis-set convergence of the relative energies is demonstrated in the SI, using basis sets up to the cc-5Z level (16,000 basis functions for the simulation cell considered). Comparing MP2 single-point energies for the structures obtained with the PBE0, PBE0-D3, and PBE approaches, we find that PBE0-D3 yields the lowest MP2 energy, making this a good starting point for subsequent optimizations. This observation is in agreement with previous work on liquid water, where it was found that PBE0-D3 closely reproduces the MP2 density.44,45 The relative energies at the MP2 and RPA level of MP2 optimized structures are displayed in Figure 4. We observe that MP2 and RPA relative energies are nearly identical and that they confirm the DFT results, with correlation factors larger than 0.995 for PBE, PBE0, and PBE0-D3 with the MP2 data. In agreement with ref 1, structures 2A1, 2B1, 9A1, and 9B2 all lie higher in energy than the experimental (2C1) structure, but we find the KSKK structure (9A2) structure is the global minimum. At all levels of theory, we find that the range of energies is ∼0.8 kJ/ mol per molecule, in good agreement with the ∼0.8 kJ/mol per molecule found in ref 9 (periodic local density functional approach) and ∼1 kJ/mol per molecule found in ref 15 (periodic semilocal density functional approach). This range is substantially less than that of >15 kJ/mol per molecule reported in ref 1. We also note non-negligible individual discrepancies between the fully periodic MP2 energy differences and those reported in ref 1 using fragment-based approaches. For example, we find that the 2C1 structure is more stable than 9B2 by 0.1 kJ/mol, whereas in ref 1 an energy difference of ∼0.8 kJ/mol is reported. The uniformity of the energy differences between structures obtained in the current work, with a range of periodic electronic structure approaches and

with two independent codes (CASTEP and CP2K), gives us confidence on the order of stability and absolute energy differences reported here. As an additional check on the sensitivity of the results to the starting geometry used, in the SI we compare the energies of structures optimized at the MP2 and RPA levels using the geometries reported in ref 1 and those obtained here and find no significant discrepancies. The density of MP2 ice XV is found to be 1.49 g/cm3, closely matching the density predicted by the PBE0-D3 results but significantly greater than the experimentally inferred value of 1.328 g/cm3 (1.476 g/cm3 for the experimentally isolated deuterated form). Looking closely at the relative energies in Figure 4, while it is clear that the experimental 2C1 structure is metastable with respect to 9A2, structures 4A2, 4C1, 4C2, and, in particular, 7C1 are quite comparable to 9A2. One possibility is that the experimental structure is a product of the slow kinetics of transformation and that annealing of the experimental sample could lead to other observed structures, such as 4C1 and 7C1. However, occupancies of the hydrogen crystallographic positions that are unambiguously associated with 2C1 increased upon cycling the sample above and below the transition temperature regime. Subsequently, further analysis, including Raman assessment,52 has shown that only the 2C1 structure can account for the diffraction and spectroscopic signatures. Because there is a clear mismatch between what is seen experimentally and what is consistently predicted by all levels of theory considered in this manuscript, we propose a possible explanation for the discrepancy. In previous studies of hydrogen order in ice, in particular, ice XI,53 and also the work of Nanda and Beran, the importance of convergence of long-range interactions has been highlighted. These electrostatic interactions are computed with an Ewald sum (see Allen and Tildesley54) or related approaches. In the common definition and implementation of the Ewald scheme, so-called tinfoil or metallic boundary conditions are employed, equivalent to assuming that the repeated unit cell is embedded in spherical cavity of a medium with infinite dielectric constant (e.g., a metal). The question here is if these metallic boundary conditions are appropriate to study the ice VI to ice XV transformation. While the influence of electrical boundary conditions has been considered in the field of classical modeling, their effects on periodic quantum mechanical calculations are comparatively rarely considered. A notable exception is recent work on ferroelectric solids, such as the work by Stengel et al.55 (and references therein), where an elegant formalism has been proposed to explore the coupling between internal energy and the applied electric field in a consistent manner. Additionally, very recently the influence of boundary conditions and depolarization fields has been discussed for ice XI.56 The hydrogen-disordered ice VI phase has a high dielectric constant, estimated by computation (after correction) to be 181 at 243 K according to ref 57, consistent with the only reported experimental value of 176 taken from ref 58. This value increases with decreasing temperature, as long as the phase remains hydrogen disordered. Ordered ice phases have low dielectric constants, around 3−5, depending on temperature. Hence, at nucleation, ice XV crystals are effectively embedded in a medium with relatively high but finite dielectric constant, and the tinfoil approximation would be reasonable. However, as a larger fraction of the ice VI transforms into an ordered phase, the effective dielectric of the medium would be reduced, and nucleation and growth could take place in a matrix of effective low dielectric constant. In line 4125



Letter

Figure 5. Relative energies in kilojoules per mole per molecule, with respect to the average, at the RI-RPA level of theory for the 18 symmetry inequivalent structures in the ice XV unit cell, including corrections for different dielectric boundary conditions.

with this model, Johari and Whalley59 reported evidence of ferroelectric ordering and a reduction in the dielectric constant of a ice VI sample with time, monitored over 252 days at a constant 117 K. However, the molecular structure of the ordered ice VI was not analyzed. To model the relative stability of ice XV nucleating within a medium of low dielectric constant, the Ewald construction requires an additional surface correction term, detailed within the works of de Leeuw et al.,60 Hummer et al.,61 and Ballenegger62 and also discussed in Boresch and Steinhauser:63 Ucorr =

2π 2 M Vf

2C1 and 9A1, 9B2 become the most stable hydrogen arrangements. Ultimately, this yields the experimental structure 2C1 essentially isoenergetic with 9A1, more stable by 0.02 kJ/ mol at the MP2 level and less stable by 0.01 kJ/mol at the RPA level, a difference that is too small to resolve by current computational procedures. However, at intermediate values of f, the 4A2, 4C1, 4C2, and 7C1 structures appear likely. As shown schematically in Figure 6, we speculate that small nuclei of

(1)

Equation 1 gives an a posteriori correction to the total energy, where M is the dipole moment of the simulation cell of volume, V, and f is (2ε + 1), when the system is embedded in a spherical cavity of a medium of dielectric constant, ε. Note that the precise expression for f depends on the shape of the cavity and will be different for platelet- or needle-shaped crystals. This term vanishes if the system has no net polarization or if the dielectric constant of the medium becomes large but will destabilize ferroelectric systems if the dielectric constant of the medium is small. To test the sensitivity of the predicted total energies of putative ice XV nuclei to the presence of a medium with finite dielectric constant, we simply employ the polarization as obtained from approximate point charges (O = −1.2, H = +0.6) on top of the MP2 geometries (reported in the SI) and add the correction (eq 1) to the computed RPA energies for f = 0, f = 190, and f = 381; we report these results in Figure 5. Note that given the large uncertainty in the dielectric constant and unknown shape of the crystallites the precise value of f is of secondary importance. However, it is immediately clear that there is a substantial influence on the relative stability of ice nuclei depending on the dielectric properties of the surrounding medium. In particular, the most striking point is that the polar structures including 9A2 are substantially destabilized by the presence of a medium with finite dielectric constant, while the structures with zero (2A1, 2B1, 2C1) and small (9A1 and 9B2) dipole moment remain unaffected. As the dielectric constant of the ice VI medium diminishes, 2A1, 2B1,

Figure 6. Schematic of the emergence of small ferroelectric KSKK nanocrystallites, whose energy depends on the dielectric properties of the surrounding mixed ice VI−XV matrix. Antiferroelectric crystals, however, can form without penalty.

ferroelectric 9A2, KSKK crystallites form within ice VI, but antiferroelectric crystallites become favored as the dielectric constant of the transforming ice VI reduces. It could be evidence for this model that experimentally only partial conversion from ice VI to ice XV is achieved of ∼78%, indicating that a proportion of the matrix is unable to transform, possibly due to the presence of ferroelectric ice XV kinetically trapped in paraelectric ice VI and antiferroelectric ice XV. Finally, other factors could also influence which hydrogen arrangement is favored. For example, DCl, which is used as a dopant to facilitate the phase transition, could influence the resultant phase, but the dopant levels used in the 4126



■

isolation of XV are rather low, approximately 1 DCl for 5500 D2O molecules. To summarize the main findings: (i) It is has been established that high-level DFT and periodic MP2 and RPA methods yield a ferroelectric Cc hydrogen-ordered structure as the lowest energy structure, in agreement with early predictions9−11 but in disagreement with more recent theoretical work1 and the experimentally observed P1̅ structure. (ii) A possible explanation for the disagreement between theory and experiment is proposed by invoking a simple model of the electrostatic environment in which ice XV crystals grow. This model suggests that complete transformation to ferroelectric ice XV is disfavored under open-circuit conditions but might be favored under closed-circuit conditions and, for example, in an applied electric field. The plausibility of (ii) could be tested directly by experiment, and indeed Jackson and Whitworth64 have reported the influence of electric fields on the formation of partially ordered ice Ih and reported evidence of the formation of a ferroelectric structure, but its detailed atomic structure was not resolved. The same mechanism could also occur in the ice Ih/XI transition and could help to explain why ice XI has never been made above ∼70% pure (in the absence of an applied field). Finally, consideration of the electrical boundary conditions might be generally important in the field of crystal structure prediction.

■

REFERENCES

(1) Nanda, K. D.; Beran, G. J. O. What Governs the Proton Ordering in Ice XV? J. Phys. Chem. Lett. 2013, 4, 3165−3169. (2) Salzmann, C. G.; Radaelli, P. G.; Slater, B.; Finney, J. L. The Polymorphism of Ice: Five Unresolved Questions. Phys. Chem. Chem. Phys. 2011, 13, 18468−18480. (3) Malkin, T. L.; Murray, B. J.; Andrey, V.; Anwar, J.; Salzmann, C. G.; Malkin, T. L.; Murray, B. J.; Brukhno, A. V.; Anwar, J.; Salzmann, C. G. Structure of Ice Crystallized from Supercooled Water. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 1041−1045. (4) Russo, J.; Romano, F.; Tanaka, H. New Metastable Form of Ice and its Role in the Homogeneous Crystallization of Water. Nat. Mater. 2014, 13, 733−739. (5) Slater, B.; Quigley, D. Crystal Nucleation: Zeroing in on Ice. Nat. Mater. 2014, 13, 670−671. (6) Quigley, D.; Alfè, D.; Slater, B. Communication: On the Stability of Ice 0, Ice i, and Ih. J. Chem. Phys. 2014, 141, 161102. (7) Sellberg, J. A.; Huang, C.; McQueen, T. a.; Loh, N. D.; Laksmono, H.; Schlesinger, D.; Sierra, R. G.; Nordlund, D.; Hampton, C. Y.; Starodub, D.; et al. Ultrafast X-ray Probing of Water Structure Below the Homogeneous Ice Nucleation Temperature. Nature 2014, 510, 381−384. (8) Salzmann, C.; Radaelli, P.; Mayer, E.; Finney, J.; Ice, X. V. A New Thermodynamically Stable Phase of Ice. Phys. Rev. Lett. 2009, 103, 105701. (9) Knight, C.; Singer, S. J. Prediction of a Phase Transition to a Hydrogen Bond Ordered Form of Ice VI. J. Phys. Chem. B 2005, 109, 21040−21046. (10) Kuo, J.-L. The Low-Temperature Proton-Ordered Phases of Ice Predicted by Ab Initio Methods. Phys. Chem. Chem. Phys. 2005, 5, 3733−3737. (11) Kuo, J.-L.; Kuhs, W. F. A First Principles Study on the Structure of Ice-VI: Static Distortion, Molecular Geometry, and Proton Ordering. J. Phys. Chem. B 2006, 110, 3697−3703. (12) Singer, S.; Kuo, J.-L.; Hirsch, T.; Knight, C.; Ojamäe, L.; Klein, M. Hydrogen-Bond Topology and the Ice VII/VIII and Ice Ih/XI Proton-Ordering Phase Transitions. Phys. Rev. Lett. 2005, 94, 135701. (13) Kuo, J.-L.; Coe, J. V.; Singer, S. J.; Band, Y. B.; Ojamäe, L. On the Use of Graph Invariants for Efficiently Generating Hydrogen Bond Topologies and Predicting Physical Properties of Water Clusters and Ice. J. Chem. Phys. 2001, 114, 2527. (14) Raza, Z.; Alfe, D.; Salzmann, C. G.; Klimeš, J.; Michaelides, A.; Slater, B. Proton Ordering in Cubic Ice and Hexagonal Ice; a Potential New Ice Phase-XIc. Phys. Chem. Chem. Phys. 2011, 13, 19788−19795. (15) Fan, X.; Bing, D.; Zhang, J.; Shen, Z.; Kuo, J.-L. Predicting the Hydrogen Bond Ordered Structures of Ice Ih, II, III, VI and Ice VII: DFT Methods with Localized Based Set. Comput. Mater. Sci. 2010, 49, S170−S175. (16) Singer, S.; Knight, C. Hydrogen Bond Ordering in Ice V and the Transition to Ice XIII. J. Chem. Phys. 2008, 129, 164513. (17) Kuo, J.-L.; Klein, M. L. Structure of Ice-VII and Ice-VIII: A Quantum Mechanical Study. J. Phys. Chem. B 2004, 108, 19634− 19639. (18) Tribello, G. A.; Slater, B. Proton Ordering Energetics in Ice Phases. Chem. Phys. Lett. 2006, 425, 246−250. (19) Singer, S. J.; Knight, C. Hydrogen-Bond Topology and Proton Ordering in Ice and Water Clusters. Adv. Chem. Phys. 2011, 147, 1−74. (20) Santra, B.; Klimeš, J.; Alfe, D.; Tkatchenko, A.; Slater, B.; Michaelides, A.; Car, R.; Scheffler, M. Hydrogen Bonds and van der Waals Forces in Ice at Ambient and High Pressures. Phys. Rev. Lett. 2011, 107. (21) Santra, B.; Klimeš, J.; Tkatchenko, A.; Alfè, D.; Slater, B.; Michaelides, A.; Car, R.; Scheffler, M. On the Accuracy of van der Waals Inclusive Density-Functional Theory Exchange-Correlation Functionals for Ice at Ambient and High Pressures. J. Chem. Phys. 2013, 139, 154702. (22) Klimeš, J.; Michaelides, A. Perspective: Advances and Challenges in Treating van der Waals Dispersion Forces in Density Functional Theory. J. Chem. Phys. 2012, 137, 120901.

ASSOCIATED CONTENT

S Supporting Information *

Detailed description of the employed methods, related computational settings, basis-set parameters, and RI-MP2 optimized geometries for all of the simulations reported. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

Letter

AUTHOR INFORMATION

Corresponding Authors

*J.V.: E-mail: [email protected]. *B.S.: E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Christoph Salzmann for important contributions to this work and a critical reading of the manuscript. B.S. thanks Gareth Tribello, Furio Cora, Julian Gale, and Andrew Walker for their contribution to previous unpublished work on ice XV. B.S. acknowledges EPSRC for funding this work via grants EP/ K038583/1 and EP/K039296/1. J.V. acknowledges financial support for this work from the European Union FP7 in the form of an ERC Starting Grant under contract no. 277910. Via our membership (B.S.) of the U.K.’s HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work made use of the facilities of HECToR and ARCHER, the U.K.’s national high-performance computing service, which is funded by the Office of Science and Technology through EPSRC’s High End Computing Programme. Calculations were also enabled by a grant from the Swiss National Supercomputer Centre (CSCS) under project ID ch5. 4127



Letter

(23) DiStasio, R. a.; Gobre, V. V.; Tkatchenko, A. Many-Body van der Waals Interactions in Molecules and Condensed Matter. J. Phys.: Condens. Matter 2014, 26, 213202. (24) Bartók, A. P.; Gillan, M. J.; Manby, F. R.; Csányi, G. MachineLearning Approach for One- and Two-Body Corrections to Density Functional Theory: Applications to Molecular and Condensed Water. Phys. Rev. B 2013, 88, 054104. (25) Gillan, M. J.; Alfè, D.; Bygrave, P. J.; Taylor, C. R.; Manby, F. R. Energy Benchmarks for Water Clusters and Ice Structures from an Embedded Many-Body Expansion. J. Chem. Phys. 2013, 139, 114101. (26) Hermann, A.; Schwerdtfeger, P. Ground-State Properties of Crystalline Ice from Periodic Hartree-Fock Calculations and a Coupled-Cluster-Based Many-Body Decomposition of the Correlation Energy. Phys. Rev. Lett. 2008, 101, 183005. (27) Booth, G. H.; Grüneis, A.; Kresse, G.; Alavi, A. Towards an Exact Description of Electronic Wavefunctions in Real Solids. Nature 2013, 493, 365−70. (28) Erba, A.; Casassa, S.; Maschio, L.; Pisani, C. DFT and LocalMP2 Periodic Study of the Structure and Stability of Two ProtonOrdered Polymorphs of Ice. J. Phys. Chem. B 2009, 113, 2347−2354. (29) He, X.; Sode, O.; Xantheas, S. S.; Hirata, S. Second-Order Many-Body Perturbation Study of Ice Ih. J. Chem. Phys. 2012, 137, 204505. (30) Macher, M.; Klimeš, J.; Franchini, C.; Kresse, G. The Random Phase Approximation Applied to Ice. J. Chem. Phys. 2014, 140, 084502. (31) Del Ben, M.; Hutter, J.; VandeVondele, J. Second-Order MøllerPlesset Perturbation Theory in the Condensed Phase: An Efficient and Massively Parallel Gaussian and Plane Waves Approach. J. Chem. Theory Comput. 2012, 8, 4177−4188. (32) Del Ben, M.; Hutter, J.; VandeVondele, J. Electron Correlation in the Condensed Phase from a Resolution of Identity Approach Based on the Gaussian and Plane Waves Scheme. J. Chem. Theory Comput. 2013, 9, 2654−2671. (33) The CP2K developers group, CP2K is freely available from: http://www.cp2k.org/ (accessed November 5, 2014). (34) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (35) Tribello, G. A.; Slater, B.; Zwijnenburg, M. A.; Bell, R. G. Isomorphism Between Ice and Silica. Phys. Chem. Chem. Phys. 2010, 12, 8597−8606. (36) Knight, C.; Singer, S. Programs to Enumerate Hydrogen Bond Topologies in Ice, And Generate Graph Invariants. Available from: https://chemistry.osu.edu/∼singer/GrEnum.html (accessed November 5, 2014). (37) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (38) Watkins, M.; Pan, D.; Wang, E. G.; Michaelides, A.; VandeVondele, J.; Slater, B. Large Variation of Vacancy Formation Energies in the Surface of Crystalline Ice. Nat. Mater. 2011, 10, 794− 798. (39) Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J.; Probert, M. I. J.; Refson, K. First Principles Methods Using CASTEP. Z. Kristollagr. 2005, 220, 567−570. (40) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0Model. J. Chem. Phys. 1999, 110, 6158−6170. (41) Schönherr, M.; Slater, B.; Hutter, J.; VandeVondele, J. Dielectric Properties of Water Ice, the Ice Ih/XI Phase Transition, and an Assessment of Density Functional Theory. J. Phys. Chem. B 2014, 118, 590−596. (42) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (43) Schmidt, J.; VandeVondele, J.; Kuo, I.-F. W.; Sebastiani, D.; Siepmann, J. I.; Hutter, J.; Mundy, C. J. Isobaric-Isothermal Molecular Dynamics Simulations Utilizing Density Functional Theory: An

Assessment of the Structure and Density of Water at Near-Ambient Conditions. J. Phys. Chem. B 2009, 113, 11959−11964. (44) Del Ben, M.; Schönherr, M.; Hutter, J.; VandeVondele, J. Bulk Liquid Water at Ambient Temperature and Pressure from MP2 Theory. J. Phys. Chem. Lett. 2013, 4, 3753−3759. (45) Del Ben, M.; Schönherr, M.; Hutter, J.; VandeVondele, J. Correction to “Bulk Liquid Water at Ambient Temperature and Pressure from MP2 Theory. J. Phys. Chem. Lett. 2014, 5, 3066−3067. (46) Lee, K.; Murray, E. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy van der Waals Density Functional. Phys. Rev. B 2010, 82, 081101. (47) Klimeš, J.; Bowler, D.; Michaelides, A. Van der Waals Density Functionals Applied to Solids. Phys. Rev. B 2011, 83, 1−13. (48) Sabatini, R.; Gorni, T.; de Gironcoli, S. Nonlocal van der Waals Density Functional Made Simple and Efficient. Phys. Rev. B 2013, 87, 041108. (49) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, And Transition Elements: Two New Functionals and Systematic Testing of four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2007, 120, 215−241. (50) Labat, F.; Pouchan, C.; Adamo, C.; Scuseria, G. E. Role of Nonlocal Exchange in Molecular Crystals: The Case of Two ProtonOrdered Phases of Ice. J. Comput. Chem. 2011, 32, 2177−2185. (51) Scuseria, G. E.; Henderson, T. M.; Sorensen, D. C. The Ground State Correlation Energy of the Random Phase Approximation from a Ring Coupled Cluster Doubles Approach. J. Chem. Phys. 2008, 129, 231101. (52) Whale, T. F.; Clark, S. J.; Finney, J. L.; Salzmann, C. G. DFTAssisted Interpretation of the Raman Spectra of Hydrogen-Ordered Ice XV. J. Raman Spectrosc. 2013, 44, 290−298. (53) Rick, S. W. Simulations of Proton Order and Disorder in Ice Ih. J. Chem. Phys. 2005, 122, 094504. (54) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: New York, 1989. (55) Stengel, M.; Spaldin, N. A.; Vanderbilt, D. Electric Displacement as the Fundamental Variable in Electronic-Structure Calculations. Nat. Phys. 2009, 5, 304−308. (56) Parkkinen, P.; Riikonen, S.; Halonen, L. Ice XI: Not That Ferroelectric. J. Phys. Chem. C 2014, DOI: 10.1021/jp510009m. (57) Aragones, J. L.; Macdowell, L. G.; Vega, C. Dielectric Constant of Ices and Water: a Lesson about Water Interactions. J. Phys. Chem. A 2011, 115, 5745−5758. (58) Johari, G. P.; Whalley, E. Dielectric Properties of Ice VI at Low Temperatures. J. Chem. Phys. 1976, 64, 4484−4489. (59) Johari, G. P.; Whalley, E. Evidence for a Very Slow Transformation in Ice VI at Low Temperatures. J. Chem. Phys. 1979, 70, 2094. (60) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Simulation of Electrostatic Systems in Periodic Boundary Conditions. I. Lattice Sums and Dielectric Constants. Proc. R. Soc. London, Ser. A 1980, 373, 27. (61) Hummer, G.; Grønbech-Jensen, N.; Neumann, M. Pressure Calculation in Polar and Charged Systems using Ewald Summation: Results for the Extended Simple Point Charge Model of Water. J. Chem. Phys. 1998, 109, 2791. (62) Ballenegger, V. Communication: On the Origin of the Surface Term in the Ewald formula. J. Chem. Phys. 2014, 140, 161102. (63) Boresch, S.; Steinhauser, O. Presumed versus Real Artifacts of the Ewald Summation Technique: The Importance of Dielectric Boundary Conditions. Ber. Bunsen-Ges. 1997, 101, 1019−1029. (64) Jackson, S. M.; Whitworth, R. W. Evidence for Ferroelectric Ordering of Ice Ih. J. Chem. Phys. 1995, 103, 7647−7648.

4128


Periodic MP2, RPA, and Boundary Condition Assessment of Hydrogen

Recommend Documents