Fifty shades of water: Benchmarking DFT functionals against

Subscriber access provided by University of Sunderland

Condensed Matter, Interfaces, and Materials

Fifty shades of water: Benchmarking DFT functionals against experimental data for ionic crystalline hydrates Getachew Kebede, Pavlin D. Mitev, Peter Broqvist, Anders Eriksson, and Kersti Hermansson J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00423 • Publication Date (Web): 31 Oct 2018 Downloaded from http://pubs.acs.org on November 6, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Fifty shades of water: Benchmarking DFT functionals against experimental data for ionic crystalline hydrates Getachew Kebede, Pavlin D. Mitev, Peter Broqvist, Anders Eriksson and Kersti Hermansson* Department of Chemistry-Ångström, Uppsala University, Box 538, SE-751 21 Uppsala, Sweden

*[email protected]

Abstract We propose that crystalline ionic hydrates constitute a valuable resource for benchmarking theoretical methods for aqueous ionic system. Many such structures are known from the experimental literature, and they represent a large variety of water-water and ion-water structural motifs. Here we have collected a dataset (CRYSTALWATER50) of 50 structurally unique "incrystal" water molecules, involved in close to 100 non-equivalent O-H∙∙∙O hydrogen bonds, and a dozen well known DFT functionals, were benchmarked with respect to their ability to describe the experimental structures and OH vibrational frequencies. We find that the PBE, RPBE-D3, and optPBE-vdW methods give the best H-bond structures and that anharmonic OH frequencies generated from B3LYP energy scans along the OH stretching coordinate with the rest of the structure taken from an optPBE-vdW-optimization (B3LYP//optPBE-vdW) outperform the other methods.

1 Environment ACS Paragon Plus

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. Introduction One of the compelling challenges when benchmarking computational methods and models is access to reliable reference data. For small molecules or clusters high-level quantum-mechanical (QM) methods such as coupled cluster (CCSD(T)) or Møller-Plesset perturbation theory (MP2, MP4) are often used as the reference. For complex condensed-matter systems, however, there exists no such golden theoretical standard, and especially none that is computationally feasible. Instead we need accurate experimental results as benchmark references. For water, a number of theoretical studies have used experimental ice and/or liquid water as benchmark systems to assess the quality of force-fields or quantum-mechanical models1-3 the ICE10 benchmark set2 for ten different ice polymorphs calculated with different DFT functionals is an example. As far as we are aware there exists no similar systematic benchmark set or assessments of water models/methods based on ionic crystalline hydrates. This is addressed in the present study. The crystalline hydrates can be seen as prototypical representatives of ionic aqueous systems, not just crystals but also ionic aqueous solutions and water/salt interfaces. Such systems all contain ion-water and hydrogen-bond interactions, and they are strongly influenced by both short- and longrange cooperative effects, arising from the ion–water • • • water and water–ion–water motifs. In short, these systems are complicated but fortunately the ionic crystalline hydrates are often structurally ordered and "well-behaved", and high-quality experimental data of many types are available. Many hundreds of ionic crystalline hydrates have been examined by single-crystal diffraction techniques and by various spectroscopies such as infrared, Raman and NMR.4-5 We propose that crystalline hydrates constitute a useful platform for benchmark studies, and in this report we use them to explore a number of common DFT functionals. A large number of exchange-correlation functionals have been developed (see for example Ref. 6) since the birth of Density Functional Theory (DFT), and a number of theories and schemes to take dispersion effects into account have been devised (see for example Refs. 7-19), not least for the bound water molecule, which is in focus in the present study. The choice of best functional to use in modelling studies is often not evident. Here we will assist in that selection, by presenting benchmark results for some selected functionals, based on results for a data-set ("CRYSTALWATER50") of 50 structurally unique "in-crystal" water molecules from 12 crystalline compounds. We will discuss issues that are important to keep in mind when selecting experimental crystalline data for benchmarking purposes. The fifty molecules examined here are involved in altogether close to 100 unique hydrogen bonds, which span an experimental R(O···O) distance range of 2.65–3.05 Å and an OH vibrational frequency range of 2775–3550 cm‒1. These broad ranges reflect that our selected crystal water collection represents a rich set of different surroundings: it includes water molecules whose nearest neighbors on the O side are monovalent, divalent, or trivalent metal cations, or water molecules, and the nearest neighbors on the H side are H-bond-accepting anions of different types, or water molecules. This diversity is illustrated in Fig. 1, which displays the nearest-neighbour surroundings for four of the fifty water molecules examined. We have calculated optimized structures, hydrogen-bond distances and in-crystal OH vibrational frequencies for the chosen compounds using 12 functionals and compare with published experimental high-quality neutron and X-ray results and infrared vibrational spectroscopy data. We will thus assess both "intermolecular" and "intramolecular" water properties. The focus of this paper is as much our analysis approach (the choice of suitable crystalline hydrates and our approach to scrutinizing the experimental data) as the outcome of our assessment of the functionals.


Page 2 of 33

Page 3 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The functionals tested belong to three families: without dispersion correction (LDA7, PBE8, RPBE9, and PBEsol10), with a posteriori dispersion corrections (PBE-D28,11, PBE-D38,12, RPBED39,12, PBEsol-D310,12) and functionals where the dispersion term is included as an integrated part of the exchange-correlation energy (vdW-DF13, vdW-DF214, optPBE-vdW15, and optB86b-vdW16). In this way we will investigate take a few steps along Jacob's ladder of DFT functionals with different exchange-correlation expressions according to Perdew and Schmidt20 and also climb a few of the lowest steps on the "stairway to heaven" as classified for dispersion-corrected functionals by Klimeš and Michaelides21. We have thus chosen the standard LDA functional with its extension to the gradient-corrected PBE density functional, and additionally two flavours of gradient-corrected functional, where one was constructed to yield an improved description of molecules and atoms (RPBE) and one of solids (PBEsol). As our crystalline hydrates are partly ionic, partly molecular, it is interesting to explore whether any of these functionals will perform better than PBE. One step up on the dispersion stairway, we will investigate how a posteriori corrections (of the Grimme type) of the dispersion interactions affect the results. In the final tests, we take one more step up and test functionals which treat the non-local correlation as an intergral part of the DFT functional together with a simultaneous improved description of the exchange energy. All in all, we believe that we have chosen a small set of functionals that allow for a systematic study of the effects of exchange and non-local correlation on the DFT results for structure and vibrational frequencies of ionic crystalline hydrates – and hopefully of other condensed-matter ionic aqueous systems as well, not yet tested here. Three recent studies of ionic crystalline hydrates in the literature are particularly relevant to mention. Thus hydrogen-bond energies were evaluated based on DFT-generated electron densities for 28 hydrogen bonds in 18 crystalline structures (four of them hydrates) in the work by Vener et al.22 , and Chaka et al.23 performed DFT calculations to study thermodynamic stabilities of a series of geologically relevant Mg-containing crystalline compounds, five of which were Mg carbonate hydrates. Although the water molecule was not in focus in either of these two studies, both of them help to demonstrate that hydrate crystals constitute excellent sources of precise information. In the context of thermal energy storage, Kiyabu et al.24 have performed high-throughput density functional theory calculations for about 160 crystalline hydrates. Although the focus in that study was rather the thermodynamics of the hydration reactions of the chosen compunds, rather than the water molecules themselves and the hydrogen bonds, the paper is nevertheless very relevant to our study as it both discusses an interesting application area of the crystalline hydrates and presents results for the optimized cell volumes with a range of density functionals.

2. Methods 2.1 Systems The experimental structures of the twelve selected crystalline hydrates25-36 are available from neutron or X-ray diffraction measurements. The space groups, the number of formula units per crystallographic unit cell (Z), and the structural references used to initiate the optimizations are all collected in Table 1. In all cases, the computational cell was chosen to be the crystallographic unit cell. The geometry optimizations were performed without imposing any symmetries or cell shape constraints. The number of structurally unique water molecules in the unit cell is given in the last column of Table 1.

2.2 Electronic calculations We performed plane-wave 3D periodic DFT calculations with 12 functionals, as listed in the introduction, using the VASP program.37-40 The valence-core electron interactions were represented 3 Environment ACS Paragon Plus


by the PAW scheme41-42 and hard pseudopotentials were used for all atoms. For H, Li, Na, Be, Mg, Al, C, N, S, O and Cl, the 1s1, 1s22s1, 2p63s1, 2s2, 2p63s2, 3s23p1, 2s22p4, 2s22p5, 3s23p4, 2s22p4 and 3s23p5 electrons, respectively, were treated as valence electrons. The k-point meshes used for the twelve hydrates were chosen as follows: a 2x3x2 Monkhorst-Pack grid for Na2CO3∙10H2O, 2x2x3 for BeSO4∙4H2O, 4x4x1 for MgSO4∙11H2O, 3x4x2 for MgCO3∙3H2O, 6x3x6 for Mg(NO3)2∙6H2O, 2x2x4 for MgSO4∙7H2O, 2x3x3 for Al(NO3)3∙9H2O, 6x3x6 for LiNO3∙3H2O, 6x6x3 for Li2SO4∙H2O, 2x2x3 for LiClO4∙3H2O, 2x3x4 for LiHCOO∙H2O and 2x2x4 for LiOH∙H2O. We used Gaussian smearing with a width of 0.1 eV in all cases. All results presented were obtained with a plane-wave kinetic energy cutoff of 1000 eV. We chose this high cut-off value after performing a convergence test for one of the crystals, Al(NO3)3∙9H2O. Our findings of the need for a high energy cut-off with the use of hard PAW pseudopotentials is consistent with the literature on water-containing systems, e.g. for the ices2. Full geometry optimizations were performed for all systems starting from the diffraction determined cell parameters and atomic positions, using neutron-diffraction structures whenever available (see the references in Table 1). It can be noted that the published cell parameters are usually obtained from powder X-ray diffraction studies, even in the experimental neutron diffraction studies, as X-ray measurement tends to be the most precise method for determination of the lattice parameters. The conjugate-gradient algorithm was used for the structure optimizations with a convergence criterion of 10–7 eV for the total energy threshold in the electronic self-consistent calculation. The cell parameters and atomic coordinates were considered converged when the forces acting on the atoms were less than 0.0021 eV/Å. Reference calculations for the isolated water molecule in a large computational box were also performed for each functional examined. 2.3 OH vibrational frequency calculations The uncoupled OH stretching frequencies of H2O(g), LiClO4∙3H2O, Li2SO4∙H2O, LiOH∙H2O and BeSO4∙4H2O were calculated using a 1-dimensional (1D) anharmonic OH vibrational model. In this model, the vibrating OH bond of water was allowed to contract and stretch while all other atomic positions remained fixed at their optimized configurations. We emphasize that using such a 1dimensional approach is not really an approximation, as we compare with isotope-isolated experimental measurements. Access to isotope-isolated experimental data has many advantages, as elaborated in the text. A reduced mass of 0.94808 amu was used for the vibrating OH bond and 19 energy points (7 and 12 energy points below and above the equilibrium OH distance for each method, re) were generated along the OH scan using a step size of 0.06 Å. The Schrödinger equation was then solved variationally for the energy levels using the discrete variable basis-set representations (DVR) approach of Light et al.43,44 Finally, the anharmonic OH stretching frequency was calculated from the energy difference between the ground and the first vibrationally excited states. This computational approach has previously been used by us to calculate OH frequencies of water45,46 and hydroxides47 in condensed H-bonded systems as well as in liquid water and ionic aqueous solutions48 and surfaces49,50.

3. Results and discussion 3.1 Structural properties Before presenting our calculated results we will make some comments about the experimental reference data selected by us. The experimental reference values for the lattice constants and atomic 4 Environment ACS Paragon Plus

Page 4 of 33

Page 5 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


positions were taken from the publications listed in Table 1. The majority of them are lowtemperature single-crystal. The agreement between our optimized cell parameters and the experimental reference data are given in Table 2. Four types of statistical averages, or "agreement indices", are given for each DFT method examined. One such agreement index is the mean absolute relative error (MARE); it is defined in the table caption, and the values are listed in the last column of the table, and illustrated in Figure 2a in terms of a dartboard plot. Hitting the origin of the dartboard means exact agreement with the experimental values (in the MARE sense). The table and the figure show that LDA and RPBE deviate more from experiment than the other methods, most of which deviate by around 1%, which is good. Among the signed indices, we find both positive and negative values which might suggest that the systematic error from the cell expansion from thermal and zero-point vibrational effects are small (see more below). Next we go on to assess the calculated hydrogen-bond distances. A number of aspects needs to be considered here, namely which H-bond distances to choose, how to deal with the temperature difference between experiment and calculation, and possibly how to deal with zero-point vibrational effects. We have chosen to monitor the R(O∙∙∙O) distances over the other common intermolecular Hbond descriptor, i.e. over the R(H∙∙∙O) distance between the donor H and acceptor O atom. This choice was made because experimental R(H∙∙∙O) distances in crystals are subject to considerable uncertainty because of the problem of locating water H-atoms using diffraction data. This is true both for X-ray diffraction (XD) and neutron diffraction (ND). It is obviously difficult to locate H atoms using XD as H is such an electron-poor atom, but also ND-determined H positions are prone to systematic errors, typically of a few hundredths of an Å (see for example the discussion in Refs. 46,51,52). The O positions determined from XD are generally well determined, both with XD and ND. Two of our experimental reference studies are from X-ray diffraction (Table 1). Next we discuss temperature effects. All the O-H···O hydrogen-bond motifs in the twelve crystals are represented in the scatter plots in Figure 3. We report the results for the experimental low-temperature and room-temperature data separately by using different colour codes (red=RT= room temperature, blue=LT=low temperature). Overall, the effect of temperature on the H-bond structure is seen to be smaller than the systematic effect of varying the functional and it has a marginal effect of the spread in the distance correlations. Having said this, one should nevertheless always try to use low-temperature structural and spectroscopic data when available. How low is low temperature? It is known that between, say, 100 and 0 K, the thermal effects on vibrational amplitudes, volumes and distances for crystalline hydrates generally tend to be quite small. To verify this we made an effort to identify single-crystal neutron diffraction studies where three different temperatures, from very low to RT, had been used. Such studies are scarce. The result of our (non-exhaustive) inventory is shown in Figure 4. For each compounds in the figure it can be noted that between the left-most point (at 20-30 K) and the second point (at 70-100 K) the distance changes are very small, less than 0.01 Å. Between ~100 K and RT, the temperature-induced change of R(O···O) is somewhat larger and depends on the particular bonding situation. The blue points in Figure 3 all refer to temperatures at, or below, 120 K except the three poins belonging to the MgCO3∙3H2O crystal which were measured at 173 K (Table 1). As for the effect of zero-point motion on the R(O···O) H-bond distances we are not aware of any experimental estimates, but Brandenburg et al.2 have presented a computational estimate of the effect for ten ice polymorphs based on dispersion-corrected HF-based free energy vs. volume scans with phonon mode contributions taken into account. Inspection of their results would suggest that the zero-point vibrations could lead to a lengthening of the experimental R(O···O) distances by up to about 0.03–0.04 Å, or expressed differently, lead to a bond contraction if one instead "back-corrects" 5 Environment ACS Paragon Plus


the experimental distances to Requilibrium(O···O) values, i.e. to the bottom of the potential energy well, as we do in calculations. Using a different approach, namely a comparison of classical and quantum-based molecular dynamics simulations (ab initio path integral MD) for a number of hydrogen-bonded clusters from the dimer to an octamer, Li et al.53 estimated the quantum effects on various properties, and found that the quantum effect led to an increase of the R(O···O) distances of 0.01 to 0.03 Å, i.e. of similar magnitude as discussed above. Clearly ZPE effects need to be considered, and extrapolating the results from these two studies to our crystalline hydrates, we note that such shifts of the point swarms in Figure 3 would be visible but small, and will hardly affect our conclusions since the shifts occurring between different functionals are (unfortunately) considerably larger than the quantum effects. Thus in our current study it appears to be meaningful to make a systematic comparison of the functionals with respect to H-bond distances also without carrying out additional ZPE-corrections for the distances. Next we discuss in some details how well (or not) the twelve functionals represent the experimental H-bond distances. Standard DFT methods, without dispersion interaction. Let us first compare the performance of local and semi-local functionals (first row in the Figure 3). It is known that LDA overestimates the H-bond strength (predicts shorter H-bond distances and larger binding energies).54,55 Also here, not surprisingly, LDA manifests its usual underestimation of the lattice constant (‒0.2 to ‒0.9 Å for the selected crystals) and H-bond distance (‒0.01 to ‒0.29 Å). The various agreement indices for the H-bond distances are given in Table 3 and the MARE values are displayed in the dartboard in Figure 2b. LDA is indeed seen to give the largest deviations from experiment, taking it out of the candidacy for good ionic crystalline hydrate studies. The GGA functionals in our collection (PBE, RPBE and PBEsol) appear to significantly improve the description of the structures of the ionic hydrates. However, RPBE and PBEsol still do not perform very well, contrary to PBE, which performs the best. As dispersion effects are not included, PBE reaches this good result by underestimating the exchange as dispersion effects are not included, and thus gets the correct answer, but for the wrong physical reasons. Therefore, without dispersion correction, or direct inclusion of non-local correlation, a standard semi-local functional with a correct exchange contribution is expected to overestimate the R(O···O) distances, just as RPBE does by a more rapid increase in the exchange with increasing density gradients. PBEsol on the other hand works in the opposite direction, and makes the results more in line with the LDA, as the exchange contribution as function of enhancement factor is damped compared to PBE. Grimme's a postriori approach. The effect of the D3 correction on PBE, RPBE and PBEsol are reported on the second row in Figure 3. The added dispersion interactions is seen to compress the cell and the R(O∙∙∙O) distances, with a noticeable effect on the agreement compared to the methods in the first row; accounting a posteriori for dispersion interaction may in practice improve or worsen the agreement between experiment and calculations. The D3 correction significantly amends the H-bond elongation tendencies of RPBE. On the other hand, as PBE already captures the R(O∙∙∙O) behaviour, adding the D3 correction places the R(O∙∙∙O) dataset "below the diagonal" in the figure and rather leads to a poorer agreement. The D2 and D3 dispersion corrections have similar influence on the results. RPBE-D3 shows the best performance in this class. Self-consistent vdW methods. Moving on to the self-consistent vdW methods, the optPBEvdW method is the best performer, giving a MARE value of 1% for the lattice constants and the data points are nicely positioned along the diagonal in the R(O∙∙∙O) plot. In our recent work regarding water on NaCl(001) and MgO(001), we showed that optPBE-vdW performs well also for the energetics of the water adsorption (this was the only quantity for which we found reference 6 Environment ACS Paragon Plus

Page 6 of 33

Page 7 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


experimental data).56 The vdW-DF functional overestimates the R(O∙∙∙O) distances. This functional has previously shown15 to predict too long H-bond distances for H-bonded complexes of the S22 data set. It is clear from our tests that the initial self-consistent vdW method vdW-DF has problems in describing hydrogen-bonded systems and that the improvements seen in the literature for other systems where dispersion is important (e.g. graphite)57 are not found here. The reason for this is believed to be due to the lack of electron density on the hydrogen. This led to the development of the vdW-DF2 density functional, but as seen, our tests with this functional did not give any significant improvement in the hydrogen bonding for the ionic hydrates. Instead, again a modified exchange description results in an improved hydrogen bond description. In our tests, optPBE-vdW stands out and gives the best result. The exchange enhancement factor in optPBE is very similar to that of PBE, so here one could again suspect that an overestimated exchange contribution corrects for the lack of the self-consistent non-local correlation effetcs. Overall, we conclude that, for the crystalline hydrates, PBE gives the best H-bond structure among the dispersion-void methods, RPBE-D3 among the a posteriori methods, and optPBE-vdW among the self-consistent vdW methods. This result is consistent with the results of Santra et al.58 and of Kiyabu et al.24, who compared cell volumes obtained with different density functionals for many ice polymorphs, and for one crystalline hydrate (MgCl2∙H2O), respectively. Next we move on to discussing an intermolecular property, ν(OH), which is one of the most used experimental measures of local structure and intermolecular interaction.

3.2 OH vibrational frequencies For crystals that contain many water molecules in the unit cell the experimental OH spectra generally consist of complex broad bands, which makes the assignment of the right frequency to the right structural feature very difficult. Such compounds are not suitable as benchmark compounds for assessment of theoretically obtained frequencies, in our opinion. For simple hydrates, however, with only few water molecules per formula unit, the experimental spectral peaks can be well resolved if the experiments have been performed at low temperature and on samples that are isotope-isolated (giving essentially uncoupled stretching frequencies). Vibrational spectra for isotope-isolated hydrates are of course much less abundant in the literature than those of the fully protonated compounds, but they do exist for a reasonably large number of crystalline hydrates. Our benchmark compounds match the requirements we have set up: (i) they are simple, (ii) high-quality lowtemperature isotope-isolated OH vibrational spectroscopic data exist for them, and (iii) high-quality experimental structure determinations have been performed. The compounds that we have chosen are LiClO4∙3H2O, Li2SO4∙H2O, BeSO4∙4H2O and LiOH∙H2O (see Tables 1 and 4). All have just one unique water molecule per unit cell. The two OH legs in the water molecule in Li2SO4∙H2O have different surroundings (Figure 1d) giving rise to two different OH frequencies, and the same is true for BeSO4∙4H2O, but for LiClO4∙3H2O and LiOH∙H2O (Figure 1c) there are symmetry elements through the water molecules which make the two OH legs in the molecule give one and the same uncoupled (isotope-isolated) frequency value. Altogether we thus have six experimental OH frequencies, plus the gas-phase frequency, in our reference set (the rightmost column of Table 4). The anharmonic OH frequencies for all these cases were calculated with six of the DFT methods used above and the resulting frequencies are compared with experiment in the first six frames in Figure 5. As usual, the diagonal means perfect agreement. The various agreement indices are given in the upper part of Table 5 and the MARE values are displayed in the dartboard in Figure 2c. 7 Environment ACS Paragon Plus


Starting with the PBE results in the upper left corner, it is clear that the calculated frequencies are seriously underestimated. This is true for all the crystals as well as for the gas-phase molecule (which is the top-most point here, and actually in all the frames). The discrepancies between the experimental and calculated frequencies is partly a consequence of discrepancies between the experimental and calculated structures although it was found in Figure 3 that, overall, the PBEoptimized structures agree well with experiment for our twelve crystals. Primarily, the underestimation of the OH vibrational frequencies is a consequence of the DFT model itself. This is seen by constructing the "Energy vs. r(OH)" potential energy scan curves underlying the vibrational calculations using the true, experimental ND-determined structure instead of the PBE-optimized structure. With this approach (not shown here), the six crystal cases still give OH frequencies that lie approximately 150-250 cm–1 below the experimental values and we can conclude that clearly it is the "Energy vs. r(OH)" scans from PBE that are too flat. The ME value for this "PBE//ND" approach, where the PBE frequency calculation was made for the ND-determined structure, is givenin the lower part of Table 5 (–199 cm–1) and is similar to the ME value from the standard PBE approach (–205 cm–1). Note that in all cases the OH distance in focus is always optimized; it is all the other atomic positions which are kept at the ND-determined values in the "PBE//ND" exercise. It is appropriate to point out that it has been observed in the literature that harmonic OH frequencies in minerals and hydrates calculated using PBE, or some similar functional, agree reasonably well with experiment. This is due to error compensation as the neglect of considering the OH anharmonicity is partially compensated by the DFT error. Moving to the right in the upper row of Figure 5 we note that the RPBE functional improves the agreement between experiment and calculations, but this is in fact a result of error cancellation. We found in the previous section that the RPBE structures are consistently too loose compared to the neutron-diffraction structures, with longer (weaker) hydrogen bonds, which one might expect would lead to a systematic upshift of the OH frequencies, thereby partly curing the large downshifts obtained with PBE. In addition, the intrinsic flattening of the "Energy vs. r(OH)" potential energy curves for RPBE is less severe than for PBE, as shown by calculating the RPBE scans based on the ND structures. The two effects work in the same direction and RPBE frequencies lie higher than PBE. For PBEsol, the "intrinsic" DFT error in the shape of the energy scan is more severe than for PBE as the scan curves are found to be considerably flatter than for PBE when evaluated using the experimental ND structures with both methods. This flatness is reflected in the large negative ME frequency for the PBEsol//ND entry in Table 5. At the same time, the optimized crystal structures are tighter, i.e. the H-bonds are shorter, for PBEsol compared to PBE or experiment; this will also suggest a lowering of the OH frequency, and indeed the total systematic error with this method becomes a very large downshift compared to experiment, as seen in the third frame of Figure 5. The too flat scan shapes generated by these three DFT methods for the crystals is also reflected in the potential energy curves for the water monomer, but to a lesser extent. For PBE the uncoupled anharmonic ν(OH) frequency for the monomer lies about 110 cm–1 below the experimental value, and this happens to be true also for PBEsol, as one can actually surmise from inspection of Figure 5 (the uppermost point is the gas-phase point in each frame). It is thus insufficient to draw conclusions from the monomer when assessing the sensitivity of the crystal frequencies to the choice of functional. As the OH stretching vibration is an intramolecular property we might expect that the dispersion should be of less importance, and that the main differences will come from the different exchange treatment, which is the main cause of the effects in the three semi-local functionals just discussed. Towards the end of this section we will test B3LYP, where part of the semi-local 8 Environment ACS Paragon Plus

Page 8 of 33

Page 9 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


exchange has been replaced with non-local Fock exchange; this is known to greatly improve the description of the chemical bonds in molecules. But first we will turn to the a posteriori dispersion correction methods, where a quick glance at Figure 5 shows that the effects of adding the Grimme D3 correction is modest. The PBE-D3 frequencies are on average very slightly more downshifted compared to experiment than the PBE frequencies. It might be tempting to expect a more substantial frequency downshift for PBE-D3 as the D3 correction was found above to give an average foreshortening of the R(O···O) distances by about 0.06 Å (Table 3), but such reasoning is not necessarily correct. For the six bound OH groups in Figure 5, the R(O···O) contractions occurring on going from PBE to PBE-D3 lie in the range 0.02 to 0.06 Å, but the corresponding OH frequency shifts are found to vary substantially, and a few even upshift in frequency. This does not contradict common wisdom: while it is well known from experiment that there exists a correlation between ν(OH) and R(O···O)59,60, the experimental correlation displays some considerable scatter around the least-squares-fitted curve. The scatter is due to the fact that R(O···O) is a reasonable but not perfect descriptor of the (full) surroundings of the vibrating OH bond; other neighbours and other structural features, such as the H-bond angle, also play major roles. Thus if we are interested in benchmarking ν(OH), it is important to study ν(OH) itself (or possibly re(OH)61), but we cannot draw firm conclusions based on a 'ν(OH) vs. R(O···O)' correlation. Now for the optPBE-vdW functional, finally, we found above that it yields good agreement with the experimental crystal structures but for the OH frequencies we should not expect good agreement as we cannot beat the systematic underestimation that plagues all the GGA-based methods. At least two different avenues can be taken to try to remedy the underestimations of OH frequencies calculated with the set of methods just discussed above. The most obvious route is to recognize that in many scientific scenarios the frequency shifts are of greater interest than the absolute frequencies themselves, as they make use of error cancellation to dampen the effects of "all sorts of" systematic errors in the absolute quantities. The gas-to-crystal frequency shift, ∆νgas=>crystal, i.e. νcrystal – νgas-phase, is indeed much used in the scientific literature to give valuable information about the magnitude of intermolecular interactions. The results for PBE, PBE-D3 and optPBE-vdW are presented in the bottom-most row in Figure 5. It is clearly seen from the dots and from the ME value given in each frame that the agreement between calculated and experimental results has improved, although the average differences are still substantial. Following another avenue, we recognize that the B3LYP hybrid functional62,63 is known to generally perform very well for molecular systems, not least for aqueous systems, and this is true also for the absolute OH frequencies. Here, optimization of the gas-phase water molecule at the B3LYP level, applying the same pseudopotentials and strict computational settings as have been used throughout this study, followed by anharmonic vibrational calculations for the uncoupled OH mode yields an OH frequency of 3704 cm–1, deviating only by 3 cm–1 from the experimental value64. Given the large computer time resources needed to optimize the systems in focus here using B3LYP, we have done the next best thing to an all-B3LYP benchmark, namely perform B3LYP frequency calculations for the optPBE-vdW-optimized structures. We use the notation B3LYP//optPBE-vdW for our approach, in line with standard notation from molecular quantum chemistry. The results for the absolute ν(OH) frequencies as well as the gas→crystal Δν(OH) frequency shifts are shown in two panels in Figure 5. The B3LYP//optPBE-vdW results are close to experiment, and about as good as the B3LYP//ND results; both are also displayed for the absolute ν(OH) frequency. In summary, the majority of DFT methods tested here give severely underestimated νcrystal and ∆νgas=>crystal values compared to the experimental reference values. optPBE-vdW gives adequate



values for ∆ν but not for νcrystal. The B3LYP//optPBE-vdW approach gives good agreement with experiment for both ∆ν and νcrystal.

4. Concluding remarks Ionic crystalline hydrates, ionic aqueous solutions and water/salt interfaces are abundant and important in nature, in biology and in a range of technical applications, where electrochemistry is just one prominent example. Consequently a very large number of modelling studies have been, and will be, performed for such systems. This highlighs the need for adequate computational models that can simultaneously handle water-water and water-ion interactions as well as cooperative effects, which are substantial in such condensed matter with polar and polarizable building blocks. Benchmark efforts are needed to assess the validity and accuracy of the interaction models for these challenging systems. We propose that ionic crystalline hydrates may constitute a suitable class of reference systems for such assessments, if chosen with care. In this study, we have selected experimental crystalline hydrate data as references for DFT calculations with twelve functionals with the purpose to assess their capabilities to reproduce two key properties of ion-water systems, namely hydrogen-bond distances and OH vibrational frequencies. Both properties are commonly used as probes of the interaction strengths, as they are sensitive to the structural details. The benchmark compounds chosen here represent a spread of different environments and span large R(O···O) and ν(OH) ranges. Our collection of hydrates contain both highly hydrated crystals and those with few water molecules per unit cell. Among the highly hydrates crystals we have MgSO4∙7H2O, MgSO4∙11H2O and Al(NO3)3∙9H2O, which are special not only in that the formula unit contains many water molecules but in that they are all structurally unique. This is something of an ultimate scenario in terms of information concentration, and it entails both scientific and practical advantages. Building a water property data-base using only a modest number of highly hydrated systems would allow us to minimize our computational effort. For the vibrational studies, on the other hand, systems with few water molecules are required so that the experiments do not yield a forest of OH peaks that are difficult to assign to the proper structural motifs. Using only vibrational data from isotope-isolated reference compounds, as is done here, removes the vibrational couplings and helps to limit the vibrational overlap but the signals from highly hydrated crystals would still be difficult to interpret. Therefore our collection consists of both water-rich and water-poor hydrates. Another issue that needs addressing in the assessment exercises is the temperature effects, namely the mismatch between the theoretical equilibrium structures obtained from the quantummechanical structure optimizations and the experimental structures which contain thermal and zeroKelvin effects. Here we made a (non-exhaustive) inventory of crystalline hydrates each of which had been studied at three temperatures in the literature: at very low T, at around 100 K and at room temperature, It was confirmed that structural data at around 100 K or below should be prioritized for benchmark structural studies. The influence of the zero-point effect on the structures was assessed by reference to two literature studies. It is worth pointing out that in our calculation of vibrational frequencies, zero-point energy effects are automatically taken into account as we use a quantumdynamical treatment. Altogether the performance of 12 density functionals was compared against experiment with respect to structures and OH frequencies for 12 periodic crystalline hydrates, containing 50 structurally unique water molecules and about one hundred hydrogen bonds. The analysis showed that the experimental ionic hydrate structures are best reproduced by the PBE, RPBE-D3, and optPBE-vdW methods. For PBE, adding the pairwise dispersion correction 10 Environment ACS Paragon Plus

Page 10 of 33

Page 11 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


term D3 compresses the cell volume, shortens R(O···O) and slightly worsens the agreement between experiment and calculations, largely a consequence of the erroneous exchange contribution in PBE. RPBE-D3 performs well, possibly for the right reasons (reasonable exchange and dispersion treatment as far as structures are concerned). The optPBE-vdW method performs the best. Our results do not contradict the results from density functional assessments of ice polymorphs in the literature, despite the fact that only very few of the fifty water molecules in our sample have only water molecules as their nearest neighbours. Santra et al.58 compared the calculated and experimental equilibrium cell volumes for seven ice polymorphs with eight functionals, two of which are the same as we have used here (PBE and optPBE-vdW). In the work of Santra et al. it turned out that, both before and after ZPE corrections, these two functionals gave the best, or among the best, volumes of all seven methods, and optPBE-vdW performed the very best. These results are consistent with our results for the ionic hydrates. Brandenburg et al.2 calculated the cell volume for 10 ice polymorphs with a large number of functionals where PBE, PBE-D3, RPBE and RPBE-D3 are in common with our study. As us, they found that RPBE severely overestimated the unit cell volumes, but that RPBE-D3, PBE and PBE-D3 gave good agreement with experiment (for PBE-D3 especially so after the ZPE back-correction of the experimental data had been undertaken in Ref. 2). Regarding the anharmonic OH frequencies, we find the B3LYP//optPBE-vdW approach, i.e. B3LYP frequency calculations on the optPBE-vdW generated structure, is an affordable and accurate approach to calculating OH frequencies for these ionic aqueous systems. It outperforms the other methods. As far as we are aware, there exist no previous systematic assessments of functionals for OH stretching vibrations in ice polymorphs or crystalline hydrates. We propose that PBE0 or other hybrid functionals with a portion of exact exchange may also work well.

Acknowledgments This work is supported by the Swedish Research Council (Vetenskapsrådet). Funding from the National Strategic e-Science program eSSENCE is greatly acknowledged. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX and NSC.



References (1) Klimeš, J.; Michaelides, A. Perspective: Advances and Challenges in Treating van der Waals Dispersion Forces in Density Functional Theory. J. Chem. Phys. 2012, 137, 120901. (2) Brandenburg, J. G.; Maas, T.; Grimme, S. Benchmarking DFT and Semiempirical Methods on Structures and Lattice Energies for Ten Ice Polymorphs. J. Chem. Phys. 2015, 142, 124104. (3) Gillan, M. J.; Alfè, D.; Michaelides, A. Perspective: How Good Is DFT for Water? J. Chem. Phys. 2016, 144, 130901. (4) Franks, F. Water: A Comprehensive Treatise. Vol.2, Water in Crystalline Hydrates, Aqueous Solutions of Simple Nonelectrolytes; Plenum Press, 1973. (5) Olovsson, I.; Jönsson, P. G. X-Ray and Neutron Diffraction Studies of Hydrogen Bonded Systems. In The Hydrogen Bond—Recent Developments in Theory and Experiments, Vol.2 , Ch. 8; Eds. P. Schuster, G. Zundel, C. Sandorfy; North-Holland, Amsterdam, 1976; pp. 393–456. (6) Marques, M. A. L.; Oliveira, M. J. T.; Burnus, T. LIBXC: A Library of Exchange and Correlation Functionals for Density Functional Theory. Comput. Phys. Commun. 2012, 183, 2272–2281. (7) Perdew, J. P.; Zunger, A. Self-Interaction Correction to Density-Functional Approximations for ManyElectron Systems. Phys. Rev. B 1981, 23, 5048–5079. (8) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (9) Hammer, B.; Hansen, L. B.; Nørskov, J. K. Improved Adsorption Energetics within Density-Functional Theory Using Revised Perdew-Burke-Ernzerhof Functionals. Phys. Rev. B 1999, 59, 7413–7421. (10) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett. 2008, 100, 136406. (11) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787–1799. (12) 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. (13) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (14) Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy van der Waals Density Functional. Phys. Rev. B 2010, 82, 081101. (15) Klimeš, J.; Bowler, D. R.; Michaelides, A. Chemical Accuracy for the van der Waals Density Functional. J. Phys. Condens. Matter Inst. Phys. J. 2010, 22, 022201. (16) Klimeš, J.; Bowler, D. R.; Michaelides, A. Van der Waals Density Functionals Applied to Solids. Phys. Rev. B 2011, 83, 195131. (17) Tkatchenko, A.; Scheffler, M. Accurate Molecular Van der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. (18) Berland, K.; Hyldgaard, P. Exchange Functional That Tests the Robustness of the Plasmon Description of the van der Waals Density Functional. Phys. Rev. B 2014, 89, 035412. (19) Bučko, T.; Lebègue, S.; Ángyán, J. G.; Hafner, J. Extending the Applicability of the TkatchenkoScheffler Dispersion Correction via Iterative Hirshfeld Partitioning. J. Chem. Phys. 2014, 141, 034114.


Page 12 of 33

Page 13 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(20) Perdew, J. P. ; Schmidt, K. In Density Functional Theory and Its Applications to Materials, Eds. V. E. Van Doren, K. Van Alsenoy and P. Geerlings, American Institute of Physics, Melville, NY, 2001. (21) Klimeš, J., Michaelides, A. Perspective: Advances and challenges in treating van der Waals dispersion forces in density functional theory. J. Chem. Phys. 2012, 137, 120901. (22) Vener, M. V.; Egorova, A. N.; Churakov, A. V.; Tsirelson, V. G. Intermolecular Hydrogen Bond Energies in Crystals Evaluated Using Electron Density Properties: DFT Computations with Periodic Boundary Conditions. J. Comput. Chem. 2012, 33, 2303–2309. (23) Chaka, A. M.; Felmy, A. R. Ab Initio Thermodynamic Model for Magnesium Carbonates and Hydrates. J. Phys. Chem. A 2014, 118, 7469–7488. (24) Kiyabu, S.; Lowe, J. S.; Ahmed, A.; Siegel, D. J. Computational Screening of Hydration Reactions for Thermal Energy Storage: New Materials and Design Rules. Chem. Mater. 2018, 30, 2006–2017. (25) Lundgren, J. O.; Liminga, R.; Tellgren, R. Neutron Diffraction Refinement of Pyroelectric Lithium Perchlorate Trihydrate. Acta Crystallogr. B 1982, 38, 15–20. (26) Tellgren, R.; Ramanujam, P. S.; Liminga, R. Hydrogen Bond Studies 78. A Neutron Diffraction Study of Lithium Formate Monohydrate, LiHCOO·H2O. Ferroelectrics 1973, 6, 191–196. (27) Hermansson, K.; Thomas, J. O. The Experimental Electron Density in Lithium Hydroxide Monohydrate. Acta Crystallogr. B 1982, 38, 2555–2563. (28) Hermansson, K.; Thomas, J. O.; Olovsson, I. Hydrogen Bond Studies. CXXXVIII. Neutron Diffraction Studies of LiNO3·3H2O at 120 and 295 K. Acta Crystallogr. B 1980, 36, 1032–1040. (29) Lundgren, J. ‐O.; Kvick, Å.; Karppinen, M.; Liminga, R.; Abrahams, S. C. Neutron Diffraction Structural Study of Pyroelectric Li2SO4·H2O at 293, 80, and 20 K. J. Chem. Phys. 1984, 80, 423–430. (30) Libowitzky, E.; Giester, G. Washing Soda (Natron), Na2CO3·10H2O, Revised: Crystal Structures at Low and Ambient Temperatures. Mineral. Petrol. 2003, 77, 177–195. (31) Kellersohn, T., Delaplane, R.G., Olovsson, I. The Synergetic Effect in Beryllium Sulfate Tetrahydratean Experimental Electron-Density Study. Acta Crystallogr. B 1984, 50, 316-326. (32) Giester, G.; Lengauer, C. L.; Rieck, B. The Crystal Structure of Nesquehonite, MgCO3 ·3H2O, from Lavrion, Greece. Mineral. Petrol. 2000, 70, 153–163. (33) Schefer, J.; Grube, M. Low Temperature Structure of Magnesium Nitrate Hexahydrate, Mg(NO3)2·6H2O: A Neutron Diffraction Study at 173 K. Mater. Res. Bull. 1995, 30, 1235–1241. (34) Fortes, A. D.; Wood, I. G.; Alfredsson, M.; Vočadlo, L.; Knight, K. S. The Thermoelastic Properties of MgSO4·7D2O (Epsomite) from Powder Neutron Diffraction and Ab Initio Calculation. Eur. J. Mineral. 2006, 18, 449–462. (35) Fortes, A. D.; Wood, I. G.; Knight, K. S. The Crystal Structure and Thermal Expansion Tensor of MgSO4·11D2O (meridianiite) Determined by Neutron Powder Diffraction. Phys. Chem. Miner. 2008, 35, 207– 221. (36) Hermansson, K. A Neutron Diffraction Determination of the Structure of Deuterated Aluminium Nitrate Nonahydrate, Al(NO3)3·9D2O. Acta Crystallogr. C 1983, 39, 925–930. (37) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558–561. (38) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the Liquid-metal-amorphousSemiconductor Transition in Germanium. Phys. Rev. B 1994, 49, 14251–14269. (39) Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15–50. (40) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186.



(41) Kresse, G.; Hafner, J. Norm-Conserving and Ultrasoft Pseudopotentials for First-Row and Transition Elements. J. Phys. Condens. Matter 1994, 6, 8245. (42) Vanderbilt, D. Soft Self-Consistent Pseudopotentials in a Generalized Eigenvalue Formalism. Phys. Rev. B 1990, 41, 7892–7895. (43) Light, J. C.; Hamilton, I. P.; Lill, J. V. Generalized Discrete Variable Approximation in Quantum Mechanics. J. Chem. Phys. 1985, 82, 1400–1409. (44) Bacic, Z.; Light, J. C. Theoretical Methods for Rovibrational States of Floppy Molecules. Annu. Rev. Phys. Chem. 1989, 40, 469–498. (45) Mitev, P. D.; Eriksson, A.; Boily, J.-F.; Hermansson, K. Vibrational Models for a Crystal with 36 Water Molecules in the Unit Cell: IR Spectra from Experiment and Calculation. Phys. Chem. Chem. Phys. 2015, 17, 10520–10531. (46) Sen, A.; Mitev, P. D.; Eriksson, A.; Hermansson, K. H-Bond and Electric Field Correlations for Water in Highly Hydrated Crystals. Int. J. Quantum Chem. 2016, 116, 67–80. (47) Hermansson, K.; Gajewski, G.; Mitev, P. D. Origin of the OH Vibrational Blue Shift in the LiOH Crystal. J. Phys. Chem. A 2008, 112, 13487–13494. (48) Hermansson, K.; Bopp, P. A.; Spångberg, D.; Pejov, L.; Bakó, I.; Mitev, P. D. The Vibrating Hydroxide Ion in Water. Chem. Phys. Lett. 2011, 514, 1–15. (49) Kebede, G. G.; Mitev, P. D.; Broqvist, P.; Kullgren, J.; Hermansson, K. Hydrogen-Bond Relations for Surface OH Species. J. Phys. Chem. C 2018, 122, 4849–4858. (50) Kebede, G. G.; Mitev, P. D.; Briels, W. J.; Hermansson, K. Red-Shifting and Blue-Shifting OH Groups on Metal Oxide Surfaces – towards a Unified Picture. Phys. Chem. Chem. Phys. 2018, 20, 12678. (51) Eriksson, A.; Berglund, B.; Tegenfeldt, J.; Lindgren, J. Corrections to the O-H Bond Lengths and H-O-H Angles of the Water Molecules in Crystalline Hydrates. Application to Ba(ClO3)2·H2O and K2C2O4·H2O. J. Mol. Struct. 1979, 52, 107–112. (52) Eriksson, A.; Hermansson, K.; Lindgren, J.; Thomas, J. O. The Use of Spectroscopic Data in the LeastSquares Refinement of Diffraction-Obtained Structure Parameters. An Example: K2C2O4·H2O. Acta Crystallogr. A 1982, 38, 138–142. (53) Li, X-Z ; Walker, B.; Michaelides, A. Quantum nature of the hydrogen bond. Proc. Natl. Acad. Sci. USA 2011, 108, 6369–6373. (54) Sim, F.; St. Amant, A.; Papai, I.; Salahub, D. R. Gaussian Density Functional Calculations on HydrogenBonded Systems. J. Am. Chem. Soc. 1992, 114, 4391–4400. (55) Laasonen, K.; Csajka, F.; Parrinello, M. Water Dimer Properties in the Gradient-Corrected Density Functional Theory. Chem. Phys. Lett. 1992, 194, 172–174. (56) Kebede, G. G.; Spångberg, D.; Mitev, P. D.; Broqvist, P.; Hermansson, K. Comparing van der Waals DFT Methods for Water on NaCl(001) and MgO(001). J. Chem. Phys. 2017, 146, 064703. (57) Chakarova-Käck, S.D.; Schröder, E.; Lundqvist, B.I. ; Langreth, D.C.Application of van der Waals Density Functional to an Extended System: Adsorption of Benzene and Naphthalene on Graphite. Phys. Rev. Lett. 2006, 96, 146107. (58) Santra, B.; Klimeš, J.; Tkatchenko, A.; Alfè, 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. (59) Berglund, B. ; Lindgren, J.; Tegenfeldt, J. On the correlation betweem deuteron quadrupole coupling constants, O-H and O-D stretching frequencies and hydrogen-bond distances in solid hydrates. J. Mol. Struct. 1978, 43, 179-191.


Page 14 of 33

Page 15 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(60) Libowitzky, E. Correlation of O-H Stretching Frequencies and O-H···O Hydrogen Bond Lengths in Minerals. Monatshefte für Chemie 1999, 130, 1047-1059. (61) Kebede, G. G.; Mitev, P. D. ; Broqvist, P.; Kullgren, J.; Hermansson, K. Hydrogen-Bond Relations for Surface OH Species. J. Phys. Chem. C 2018, 122, 4849−4858. (62) Becke, A. D. Density‐functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. (63) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785–789. (64) Benedict, W. S.; Gailar, N.; Plyler, E. K. Rotation‐Vibration Spectra of Deuterated Water Vapor. J. Chem. Phys. 1956, 24, 1139–1165. (65) Berglund, B.; Lindgren, J.; Tegenfeldt, J. O-H and O-D Stretching Vibrations in Isotopically Dilute HDO Molecules in Some Solid Hydrates. J. Mol. Struct. 1978, 43, 169–177. (66) Hayward, H. P.; Schiffer, J. Intramolecular Forces in the Water Molecule in Lithium Sulfate Monohydrate. J. Chem. Phys. 1975, 62, 1473–1479. (67) Kvick, Å.; Liminga, R.; Abrahams, S.C. Neutron diffraction study of Ba(NO2)2·H2O at 298, 102, and 20 K. J. Chem. Phys. 1982, 76, 5508-5514. (68) Küppers, H.; Takusagawa, F.; Koetzle, T.F. Neutron diffraction study of lithium hydrogen phthalate monohydrate: A material with two very short intramolecular O···H···O hydrogen bonds. J. Chem. Phys. 1985, 82, 5636-5647.

(69) Lundgren, J.-O.; Kvick, Å.; Liminga, R.; Abrahams, S.C. Neutron diffraction structural study of pyroelectric Sr(NO2)2·H2O at 294, 100, and 20 K. J. Chem. Phys. 1985, 83, 2426-2434.



Page 16 of 33

Table 1. The crystalline ionic hydrate systems studied in this work together with their space groups, number of formula units per crystallographic unit cell (Z), diffraction technique [X-ray diffraction (XD) or neutron diffraction (ND)], samples (single crystal or powder), and temperature, and the experimental references from which our starting structures were taken.

Space group

Hydrate

Z

Expt. reference

Technique, sample, Temp

No. of structurally unique H2O

LiClO4∙3H2O

P63mc

2

Lundgren et al. (1982) Ref. 25

ND, single x-l, 294 K

1

LiHCOO∙H2O

Pna21

4

Tellgren et al. (1973) Ref. 26


1

LiOH∙H2O

C2/m

4

Hermansson & Thomas (1982) Ref. 27


1

LiNO3∙3H2O

Cmcm

4

Hermansson et al. (1980) Ref. 28


2

Li2SO4∙H2O

P21

2

Lundgren et al. (1984) Ref. 29


1

Na2CO3∙10H2O

Cc

4

Libowitzky & Giester (2003) Ref. 30

XD, single x-l, 110 K

10

BeSO4∙4H2O

I4c2

4

Kellersohn (1994) Ref. 31


1

MgCO3∙3H2O

P21/n

4

Giester et al. (2000) Ref. 32

XD, single x-l, 300 K

3

P21/c

2

Schefer et al. (1995) Ref. 33

ND, single x-l,173 K

3

P212121

4

Fortes et al. (2006) Ref..34

ND, powder, 2 K

7

P1

4

Fortes et al. (2008) Ref. 35

ND, powder, 4.2 K

11

P21/c

4

Hermansson (1983) Ref. 36


Mg(NO3)2∙6H2O MgSO4∙7H2O

a

MgSO4∙11H2O

a

Al(NO3)3∙9H2O

a

Total a The diffraction measurements were made for deuterated crystals.


9 50

Page 17 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2. Errors in the optimized cell parameters obtained with the different functionals tested in this study; the errors were calculated with respect to the experimental references structures listed in Table 1. The mean error (ME), mean absolute error (MAE), mean relative error (MRE) and mean absolute relative error (MARE) were calculated using the experimental lattice constants as references (with N=36). The errors listed are defined as 1

1

𝑁 𝑁 𝑀𝐸 = 𝑁∑𝑖 = 1𝐷𝐹𝑇 ― 𝐸𝑥𝑝 , 𝑀𝐴𝐸 = 𝑁∑𝑖 = 1|𝐷𝐹𝑇 ― 𝐸𝑥𝑝| ,

𝑀𝑅𝐸 =

100 𝑁 𝐷𝐹𝑇 ― 𝐸𝑥𝑝 ∑ , 𝑁 𝑖=1 𝐸𝑥𝑝

Statistical errors for the cell parameters ME (Å)

MAE (Å)

MRE (%)

MARE (%)

LDA

‒0.32

0.32

‒3.7

3.7

PBE

+0.10

0.12

+1.3

1.5

RPBE

+0.35

0.35

+4.2

4.2

PBEsol

‒0.07

0.12

‒0.8

1.4

PBE-D2

‒0.06

0.10

‒0.8

1.2

PBE-D3

‒0.06

0.09

‒0.7

1.0

RPBE-D3

+0.07

0.09

+0.9

1.1

PBEsol-D3

‒0.12

0.20

‒2.3

2.4

vdW-DF

+0.17

0.17

+2.0

2.0

vdW-DF2

+0.98

0.11

+1.1

1.2

optPBE-vdW

+0.04

0.07

+0.4

0.8

optB86b-vdW

‒0.06

0.08

‒0.7

1.0

Method


𝑀𝐴𝑅𝐸 =

|

|

100 𝑁 𝐷𝐹𝑇 ― 𝐸𝑥𝑝 ∑ 𝑁 𝑖=1 𝐸𝑥𝑝


Table 3. Errors in the calculated R(O∙∙∙O) distances for the optimized structures compared to the experimental references structures listed in Table 1. The definitions of the errors are given in Table 2, but here N=70 for the comparisons with low-temperature (LT) experimental data and N=28 for the comparisons with room-temperature (RT) data. Statistical errors for the R(O∙∙∙O) distances: LT Method

ME (Å)

MAE (Å)

MRE (%)

MARE (%)

LDA

-0.13

0.13

-4.7

4.7

PBE

+0.02

0.03

+0.6

1.6

RPBE

+0.05

0.09

+1.8

4.7

PBEsol

-0.07

0.05

-2.3

2.7

PBE-D2

-0.03

0.03

-1.1

1.7

PBE-D3

-0.04

0.04

-1.5

1.9

RPBE-D3

+0.02

0.02

+0.6

1.0

PBEsol-D3

-0.09

0.06

-3.4

3.4

vdW-DF

+0.06

0.04

+2.2

2.2

vdW-DF2

+0.04

0.03

+1.3

1.6

optPBE-vdW

+0.01

0.02

+0.4

0.9

optB86b-vdW

-0.03

0.03

-1.0

1.3

Statistical errors for the R(O∙∙∙O) distances: RT Method

ME (Å)

MAE (Å)

MRE (%)

MARE (%)

LDA

-0.14

0.14

-5.0

5.0

PBE

-0.02

0.03

-0.6

1.0

RPBE

+0.07

0.08

+2.6

2.7

PBEsol

-0.08

0.08

-3.0

3.0

PBE-D2

-0.06

0.06

-2.1

2.1

PBE-D3

-0.05

0.05

-1.9

1.9

RPBE-D3

-0.01

0.02

-0.3

0.6

PBEsol-D3

-0.11

0.11

-3.8

3.8

vdW-DF

+0.04

0.04

+1.6

1.6

vdW-DF2

+0.02

0.03

+0.9

1.0

optPBE-vdW

-0.01

0.02

-0.3

0.6

optB86b-vdW

-0.05

0.05

-1.7

1.8


Page 18 of 33

Page 19 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 4. The crystalline ionic hydrate systems with adhering experimental OH stretching frequencies used as references for benchmarking in this study.

Hydrate

Expt. reference

Technique, isotopic ratio, Temp.

Isotope-isolated OH frequencies (cm–1)

HDO(g)

Benedict et al. (1956) Ref. 64

Microwave spectr.

3707

LiClO4∙3H2O

Berglund et al. (1978) Ref. 65

IR, 5.0 % H2O, 90 K

3556

LiOH∙H2O


IR, 5.0 % H2O, 90 K

2775

Li2SO4∙H2O

Hayward & Schiffer (1975) Ref. 66

IR, 5.4 % H2O, 20 K

3472, 3477

BeSO4∙4H2O


IR, 5.0 % H2O, 90 K

2906, 3145



Table 5. Errors in the calculated anharmonic OH frequencies, ν(OH), calculated for the selection of crystalline hydrate structures listed in Table 4. The definitions of the errors are given in Table 2, but here N=7. B3LYP//optPBE-vdW means that the potential energy scan for the frequency calculations was performed at the B3LYP level for the optPBE-vdW optimized structure. The //ND notation means that the potential energy scans were performed for the neutron-diffraction-determined structures taken from the literature (structural references in Table 1).

Statistical errors for the ν(OH), fully optimized structure ME (cm–1)

MAE (cm–1)

MRE (%)

MARE (%)

PBE

‒205

205

‒6.4

6.4

RPBE

‒145

145

‒4.6

4.6

PBEsol

‒280

280

–8.7

8.7

PBE-D3

‒219

219

‒6.9

6.9

RPBE-D3

‒176

176

‒5.6

5.6

optPBE-vdW

‒189

189

‒5.9

5.9

B3LYP//optPBE-vdW

‒38

38

‒1.2

1.2

Method

Statistical errors for the ν(OH), based on ND structiure ME (cm–1)

MAE (cm–1)

MRE (%)

MARE (%)

PBE//ND

-199

199

-6,3

6,3

RPBE//ND

-167

167

-5,3

5,3

PBEsol//ND

-243

243

-7,7

7,7

PBE-D3//ND

-207

207

-6,5

6,5

RPBE-D3//ND

-182

182

-5,8

5,8

optPBE-vdW//ND

-201

201

-6,3

6,3

B3LYP//ND

-51

51

-1,6

1,6

Method


Page 20 of 33

Page 21 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure captions Figure 1. Coordination figures around 4 of the 50 water molecules investigated. Figure 2. Dartboard plots of the Mean Absolute Relative Error (MARE) values, as defined in the caption of Table 2 and listed in the right-most columns in Tables 2 and 3. (a) Cell parameters, (b) R(O∙∙∙O), and (c) ν(OH). Figure 3. Calculated vs. experimental R(O∙∙∙O) distances. See Table 1 for the experimental reference publications. For some of the frames, a few (≤ 4) calculated points lie above the frame; all belonging to one particular water molecule which forms bifurcated hydrogen bonds, namely in LiNO3·3H2O. A few odd points which are seen to deviate very much from the diagonal line for PBEsol, PBE-D3 and PBEsol-D3 also refer to this water molecule. Figure 4. Experimental literature data for the temperature dependence of R(O∙∙∙O) distances from five crystalline hydrates whose structures have been determined from neutron diffraction data at three different temperatures; one at very low temperature, one around 100 K and one around room temperature. The references of the diffraction studies are Ba(NO2)2·H2O (Ref. 67), Li2SO4·H2O (Ref. 29), LiC8H5O4·H2O (Ref. 68), Sr(NO2)2·H2O (Ref. 69), and BeSO4·4H2O (Reg. 31). Figure 5. Calculated vs. experimental absolute OH frequencies (ν(OH), two upper rows) and gasto-crystal frequency shifts (Δν(OH), third row). See Table 4 for the experimental reference values. The wavenumber written in each frame in the upper two rows is the Mean error (ME), taken from Table 5. The ME values for the Δν(OH) in the bottom row also make use of the gas-phase frequencies.



(a) W2 in Al(NO3)3•9H2O(s)

(c) W in LiOH•H2O(s)

Page 22 of 33

(b) W7 in Al(NO3)3•9H2O(s)

(d) W in Li2SO4•H2O(s)

Figure 1 22 Environment ACS Paragon Plus

Page 23 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a) Cell volumes

(b) R(O•••O)

(c) ν(OH)

Figure 2



Page 24 of 33

Figure 3


Page 25 of 33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4



–205 cm–1

–145 cm–1

–280 cm–1

Page 26 of 33

–219 cm–1

2 points

–176 cm–1

–189 cm–1

–110 cm–1

–124 cm–1

–38 cm–1

–54 cm–1

–51 cm–1

–20 cm–1

Figure 5


Page 27 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


“For Table of Contents Only”

50 unique water molecules, 100 unique H-bonds



Figure 1. Coordination figures around 4 of the 50 water molecules investigated. 191x182mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 28 of 33

Page 29 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Dartboard plots of the Mean Absolute Relative Error (MARE) values, as defined in the caption of Table 2 and listed in the right-most columns in Tables 2 and 3. (a) Cell parameters, (b) R(O∙∙∙O), and (c) ν(OH). 771x250mm (300 x 300 DPI)



Figure 3. Calculated vs. experimental R(O∙∙∙O) distances. See Table 1 for the experimental reference publications. For some of the frames, a few (≤ 4) calculated points lie above the frame; all belonging to one particular water molecule which forms bifurcated hydrogen bonds, namely in LiNO3•3H2O. A few odd points which are seen to deviate very much from the diagonal line for PBEsol, PBE-D3 and PBEsol-D3 also refer to this water molecule. 193x161mm (300 x 300 DPI)


Page 30 of 33

Page 31 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


392x392mm (72 x 72 DPI)



Figure 5. Calculated vs. experimental absolute OH frequencies (ν(OH), two upper rows) and gas-to-crystal frequency shifts (Δν(OH), third row). See Table 4 for the experimental reference values. The wavenumber written in each frame in the upper two rows is the Mean error (ME), taken from Table 5. The ME values for the Δν(OH) in the bottom row also make use of the gas-phase frequencies. 180x160mm (300 x 300 DPI)


Page 32 of 33

Page 33 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


“For Table of Contents Only” 89x35mm (300 x 300 DPI)


Fifty shades of water: Benchmarking DFT functionals against

Recommend Documents